% Consider Example 1 in Section 5.5 on page 279 of the book
% This is a physical system consisting of 2 masses and two springs
% that we considered in class.
% We have m1 = 2, m2 = 1, k1 = 4, and k2 = 2.
% Initial conditions are x(0) = 3, x'(0) = 0, y(0) = 3, and y'(0) = 0.
% We want to find expressions for x(t) and y(t)
% Step 1: Write down a system of first order differential equations
%  that gives the state space representation of the physical system.
% Let x1 = x, x2 = x', x3 = y, and x4 = y'
% We got the following representation (see class notes)
% X' = AX
% where A = [0 1 0 0; -(k1+k2)/m1 0 k2/m1 0; 0 0 0 1; k2/m2 0 -k2/m2 0]
% Substituing system values we obtain
% A = [0 1 0 0; -3 0 1 0; 0 0 0 1; 2 0 -2 0];
% Step 2: Compute the eigenvalues and eigenvectors of A in MATLAB

A = [0 1 0 0; -3 0 1 0; 0 0 0 1; 2 0 -2 0];
[V,D] = eig(A)

% MATLAB returns the following result

V =

        0 - 0.3162i        0 + 0.3162i   0.3162             0.3162          
   0.6325             0.6325                  0 + 0.3162i        0 - 0.3162i
        0 + 0.3162i        0 - 0.3162i   0.6325             0.6325          
  -0.6325            -0.6325                  0 + 0.6325i        0 - 0.6325i


D =

        0 + 2.0000i        0                  0                  0          
        0                  0 - 2.0000i        0                  0          
        0                  0                  0 + 1.0000i        0          
        0                  0                  0                  0 -1.0000i

% So our four eigenvalues are purely imaginary. This is to be expected
% since there is no physical damping in the system
% Our eigenvalues are r1 = 2i, r2 = -2i, r3 =i, and r4 = -i
% The associated eigenvectors v1,v2,v3, and v4 are given by the
% corresponding columns of the V matrix.
%  For example v1 = [-0.3162i; 0.6325; 0.3162i; -0.6325] etc.

% Write down the general solution X(t) for the system using the
% eiganvalue and eigenvector information 
% See the discussion on complex eigenvalues in Section 9.6 on page 546
% of the book.
% For the complex conjugate eigenvalues 2i and -2i we have
% alpha = 0 and beta = 2
% Also the corresponding eigenvectors are a+ib, a-ib where
% a = [0; 0.6325; 0; -0.6325] and b = [-0.3162; 0; 0.3162; 0]
%  Use formulas (6) and (7) to find X1(t) and X2(t).
% Repeat the discussion for the complex conjugate eigenvalues i and -i
% to obtain X3(t) and X4(t)
% Now X(t) = c1*X1(t) + c2*X2(t) + c3*X3(t) + c4*X4(t)
% I got X(t) = B(t)*c where
%  X(t) = [x1(t); x2(t); x3(t); x4(t)]
% B(t) = [0.3162*sin(2t) -0.3162*cos(2t) 0.3162*cos(t) 0.3162*sin(t);
%           0.6325*cos(2t) 0.6325*sin(2t)  -0.3162*sin(t) 0.3162*cos(t);            
%         -0.3162*sin(2t)  0.3162*cos(2t)  0.6325*cos(t) 0.6325*sin(t);
%        -0.6325*cos(2t) -0.6325*sin(2t) -0.6325*sin(t) 0.6325*cos(t)]
% and c = [c1; c2; c3; c4]
% We have X(0) = [3; 0; 3; 0]
% Now X(0) = B(0)*c gives
% d = X(0) = [3; 0; 3; 0]
% B = B(0) = [0 -0.3162 0.3162 0; 0.6325 0 0 0.3162; 0 0.3162 0.6325 0;
%                  -0.6325 0 0 0.6325]
% Let us solve B*c = d in MATLAB

B = [0 -0.3162 0.3162 0; 0.6325 0 0 0.3162; 0 0.3162 0.6325 0; 
-0.6325 0 0 0.6325];
d = [3; 0; 3; 0];
c = B\d

% MATLAB returns the following answer

c =

         0
   -3.1632
    6.3244
         0

% Therefore X(t) = B(t)*c with this value of c
% Carrying out the matrix vector multiplication gives
% x1(t) = x(t) = 2*cos(t) + cos(2t) and
% x3(t) = y(t) = 4*cos(t) - cos(2t)
% This agrees with the answer given at the end of page 281 of the book
% The plots of x(t) and y(t) are given in Figure 5.22 on page 282 of
% the book.
% Enjoy!