The solution to the vector differential equation $\mathbf y'(x)=A\mathbf y(x)$ is, not surprisingly, $e^{xA}\mathbf C$, where $\mathbf C$ is a vector of constants determined by boundary conditions. The exponential of a matrix is defined via a power series, but in practice one doesn’t use that to compute it.
If $A$ is diagonalizable, it can be decomposed as $B\Lambda B^{-1}$, where $B=\begin{bmatrix}\mathbf b_1,\cdots,\mathbf b_n\end{bmatrix}$ with eigenvectors of $A$ as its columns and $\Lambda=\operatorname{diag}(\lambda_1,\dots,\lambda_n)$ with $\lambda_j$ the eigenvalues corresponding to $\mathbf b_j$. Just as $A^k=B\Lambda^kB^{-1}=B\operatorname{diag}(\lambda_1^k,\dots,\lambda_n^k)B^{-1}$, so, too $e^{xA}=Be^{x\Lambda}B^{-1}=B\operatorname{diag}(e^{\lambda_1x},\dots,e^{\lambda_nx})B^{-1}$. So, the solution to the equation is $$e^{xA}\mathbf C=Be^{x\Lambda}B^{-1}\mathbf C.$$ Since $B$ has full rank, $B^{-1}\mathbf C$ is also a vector of arbitrary constants, so we can expand the above expression as $$C_1e^{\lambda_1x}\mathbf b_1+\cdots+C_ne^{\lambda_nx}\mathbf b^n.\tag1$$
In this problem, $$A=\begin{bmatrix}-2&1&-2\\1&-2&2\\3&-3&5\end{bmatrix}.$$ Its eigenvalues can be found to be $3$, $-1$ and $-1$. We compute the kernel of $3I-A$ via row-reduction: $$\begin{bmatrix}5&-1&2\\-1&5&-2\\-3&3&-2\end{bmatrix}\to\begin{bmatrix}1&0&\frac13\\0&1&-\frac13\\0&0&0\end{bmatrix}$$ so an eigenvector of $3$ is $(-1,1,3)^T$. Moving to $A+I$, $$\begin{bmatrix}-1&1&-2\\1&-1&2\\3&-3&6\end{bmatrix}\to\begin{bmatrix}1&-1&2\\0&0&0\\0&0&0\end{bmatrix}$$ from which we get $(1,1,0)^T$ and $(-2,0,1)^T$ as linearly independent eigenvectors. Plugging these values into (1), the solution to the equation is therefore $$C_1e^{3x}\begin{bmatrix}-1\\1\\3\end{bmatrix}+C_2e^{-x}\begin{bmatrix}1\\1\\0\end{bmatrix}+C_3e^{-x}\begin{bmatrix}-2\\0\\1\end{bmatrix}=\begin{bmatrix}-C_1e^{3x}+(C_2-2C_3)e^{-x}\\C_1e^{3x}+C_2e^{-x}\\3C_1e^{3x}+C_3e^{-x}\end{bmatrix}.$$
From the image you can read of a period of about $6$ which confirms the label of a frequency $ω_0=1$. Thus the basis is a harmonic oscillator $y''+(y-1)=0$. In addition there is a friction term that leads to the represented exponential decay, where one can guess that the friction coefficient is $2\zeta$ (as per the normalized form in the answer of Mostafa Ayaz). In total this makes an equation
$$y''+2\zeta y'+y=x.$$
As per usual in such situations of exploring system responses it is assumed that all functions are zero for $t<0$, so that the initial conditions are $y(0)=0$ and $y'(0)=0$. The unit step function is $x(t)=1$ for $t\ge 0$, and of course $x(t)=0$ for $t<0$. Then solve
$$y''+2\zeta y'+y=1$$
for $t\ge 0$.
t = np.linspace(0,30,1501)
for zeta in [0.1,0.2,0.5,1,2,4]:
sol = odeint(lambda y,t:[y[1], 1-y[0]-2*zeta*y[1]],[0,0], t)
plt.plot(t,sol[:,0], label="$\zeta=$%.2f"%zeta);
plt.grid(); plt.legend(); plt.show()
to get
Best Answer
Just consider parametric curves of the form $$ (y_1(t),y_2(t)), t \in \mathbb{R} $$
For each set of initial conditions you will obtain a stream line. The picture below was obtained with mathematica using the StreamPlot[] command. In red one of the curves I mentioned before.