Two coupled nonlinear mechanical oscillators are represented by the set of differential equations:
$$\ddot{x} + \xi_1 \dot{x} + x - \gamma_1 z-\gamma_3z^3 - \xi_2\dot{z}=f(t)$$
$$\mu \ddot{z}+\gamma_1 z + \gamma_3 z^3 + \xi_2\dot{z}=-\mu \ddot{x}$$
both $x$ and $z$ are functions of time ($x(t)$, $z(t)$) and the overdot is the time derivative.
I am seeking a method to obtain only one ODE (4th order), which can be written in terms of $x$ or $z$, using the elimination method.
I have tried the operator method writing the two equations as
$$ \left(D^2 + \xi_1 D + 1\right)x - \left(\gamma_1 + \gamma_3 z^2 +\xi_2 D \right) z = f(t)$$ $$ \left(\mu D^2 + \gamma_1 + \gamma_3 z^2 + \xi_2 D\right)z + \mu D^2 x = 0$$
where $D = \frac{d}{dt}$
How to proceed from here? (Can this be done using computer algebra software ?)
Add the two equations together.
$$(1+\mu)x''+\xi_1x'+x = -\mu z''+ f(t)$$
Notice, that $(1+\mu)[x''+\frac{\xi_1}{1+\mu}x'+{1\over{1+\mu}}x]= (1+\mu)(D+\alpha)(D+\beta)x$, where $\alpha + \beta = \frac{\xi_1}{1+\mu}$ and $\alpha\beta={1\over{1+\mu}}$.
Now $(D+c)g=e^{-ct}\frac{d}{dt}(e^{ct}g)$. So, expanding the $lhs$ above, we get
$$(1+\mu)e^{-\alpha t}\frac{d}{dt}\left(e^{\alpha t} \left[e^{-\beta t}\frac{d}{dt}(e^{\beta t} x) \right] \right)=(1+\mu)e^{-\alpha t}\frac{d}{dt}\left( e^{(\alpha - \beta )t} \frac{d}{dt}(e^{\beta t} x) \right)$$
Substituting this back into the original equation we see $$(1+\mu)e^{-\alpha t}\frac{d}{dt}\left( e^{(\alpha - \beta) t} \frac{d}{dt}(e^{\beta t} x) \right) = -\mu z'' + f(t)$$
We solve for $x$ as follows $$e^{(\alpha - \beta) t} \frac{d}{dt}(e^{\beta t} x)= \int{}-\frac{\mu}{1+\mu}e^{\alpha t}z'' + {1\over{1+\mu}}e^{\alpha t}f dt.$$
Using integration by parts, $$\begin{align} \int{}-\frac{\mu}{1+\mu}e^{\alpha t}z'' dt &= -\frac{\mu}{1+\mu}\left(z'e^{\alpha t}-\int \alpha e^{\alpha t} z'dt\right) \\ &=-\frac{\mu}{1+\mu}\left(z'e^{\alpha t}-\left[\alpha e^{\alpha t}z-\int \alpha^2 e^{\alpha t} z dt\right]\right) \end{align}$$ Simplifying, $$\frac{d}{dt}(e^{\beta t} x) = -\frac{\mu}{1+\mu}e^{\beta t}\left(z'- \alpha z + \int \alpha^2 e^{\alpha (\tau - t)} z d\tau - {1\over\mu} \int e^{\alpha (\tau - t)} f d\tau\right)$$ In similar fashion, we can do the final integration and find $x$ in terms of $z$ and $f$.
Because I do not want to bore you with the derivation, let me show you how you could get a similar result with Laplace's transform. Let $X$,$Z$,$F$ be transforms of $x$,$z$,$f$ respectively. Then, the first equation becomes $$\begin{align} \left[(1+\mu)s^2+\xi_1s+1\right]X &= -\mu s^2Z+ F\\ X &= \frac{-\mu s^2}{(1+\mu)s^2+\xi_1s+1}Z + \frac{1}{(1+\mu)s^2+\xi_1s+1}F \\ X &= \frac{-\mu}{1+\mu} \frac{s^2 Z}{(s+\alpha)(s+\beta)}+\frac{1}{1+\mu} \frac{F}{(s+\alpha)(s+\beta)} \\ X &= \frac{-\mu}{1+\mu}\left(1+\frac{1}{\beta-\alpha}\left[\frac{\alpha^2}{s+\alpha}-\frac{\beta^2}{s+\beta}\right]\right)Z + \\ &\frac{1}{1+\mu} \frac{1}{\beta-\alpha}\left({1\over{s+\alpha}}-{1\over{s+\beta}}\right)F \end{align}$$
You can then apply the inverse transform of terms of the form ${1\over{s+c}}Y$, which is $e^{ct}\ast y(t)$, with appropriate assumptions about the support of the signals, namely that they don't extend from 0 to $-\infty$.