Our system is a two DOF system, $\newcommand\bs\boldsymbol$ its structural matrices are

$$ \boldsymbol{M} = m\,\begin{bmatrix}2&0\\0&1\end{bmatrix},\qquad \boldsymbol{K} = k\,\begin{bmatrix}3&-2\\-2&2\end{bmatrix}. $$

The system is harmonically loaded,

$$ \boldsymbol{p}(t)=p_o\begin{Bmatrix}0\\1\end{Bmatrix}\sin\omega t=p_o\,\bs {r}\sin\omega t, $$

where we have introduced an adimensional load shape vector.

A particular integral can be $\bs{x}_{ss}=\bs{\xi}\sin\omega t$, substituting in the equation of motion and simplifying the time dependency

$$(k \bs K -\omega^2m \bs M ) \, \bs\xi= p_o \, \bs r.$$

Introducing the unit frequency, defined in terms of unit stiffness and unit mass, $\omega_o = \sqrt{\frac km}$, with $\omega = 2 \omega_o$, dividing both members by $k$ we have

$$(\bs K - 4 \bs M) \, \bs\xi = \frac{p_o}k \,\bs r = \Delta \,\bs r.$$

Substituting the numerical values, solving for $\bs\xi$ and substituting in $\bs x_s$

$$\bs{x}_s(t) = \Delta \frac16\begin{Bmatrix}+2\\-5\end{Bmatrix}\sin\omega t.$$

If our systems starts from rest conditions, the steady state solution has a non zero initial velocity, so we have to superpose a homogeneous solution to get the respect of all initial conditions.

To find the homogeneous solution we use separation of variables, $\bs x = \bs\psi\sin\omega t$, substituting $x(t)$ in the equation of free vibrations and simplifying the time dependency we have the following homogeneous equation

$$(\bs K -\omega^2\bs M)\bs\psi=\bs0.$$

The non trivial solutions can be found using the library function eigh, that computes an 1-D array of eigenvalues and a 2-D array of (mass-normalized) eigenvectors.

Note that both evals and evecs are adimensional, the dimensional eigenvalues can be obtained multiplying evals by $\omega_o^2$.

While we are at it, we compute also the adimensional frequencies, $\omega_i/\omega_o$, and the inverses of the adimensional frequencies.

The first eigenvector has both components negative, I don't like that...

Now, the load vector divided $p_o$ and the excitation frequency, $\omega_f/\omega_o$:

First, we compute $\bs\xi/\Delta$, then its components in the modal coordinates...

The steady state response, in modal coordinates, is

$$x_s(t) = \Delta\,(\bs\psi_1 q_{s,1} + \bs\psi_2 q_{s,2})\sin\omega t$$

and for initial rest conditions we have

$$q_i(t) = \Delta \, \left(\sin\omega t - \frac{\omega}{\omega_i}\sin\omega_i t\right) q_{s,i}.$$

In code:

It's high time to represent our results, here the response is plotted in terms of adimensional modal coordinates vs adimensional time.

It is worth noting that the second frequency is rather close to the frequency of the harmonmic loading, and this is clearly reflected in the beating behaviour of the second modal coordinate.

Now, we define quite naively the response in natural coordinates and proceed to plotting in a manner similar to our previous graph.