% Using eigenvalues and vectors to solve a system of ODEs

A  =  [0 -1; 1 0];    % matrix that defines the linear ODEs
x0 =  [1 0]';         % initial condition
t=3*pi/2;             % set the time at which you want to evaluate x

[V,D] = eig(A);       % find eigenvalues and vectors of A

 
% The solution is given by the  command:

x_of_t = V  *  diag( exp( diag(D*t) ) ) * ( V \ x0 )

% ( V \ x0 )  solves for the arbitrary solution constants in terms of ICs.

% diag( exp( diag(D*t) ) )  gives   exp(lambda*t) on the diagonals
% and zero elsewhere. 

% The whole command thus gives a linear combination of the eigenvectors
% of A with each eigenvector multiplied by the appropriate constant 
% to satisfy the ICs and by e raised to the power of (the associated 
% eigenvalue lambda)*t  Note that the solution is real but several of 
% the intermediate calculations are complex.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% NEWS FLASH (12/4/96): The command mistakenly used in both Ruina
% and Terrell's lectures as well as in Nagy's solutions to 
% homework 8 (which was mistitled 'Solution to HW#7') was:

x_of_t = V  *   exp( D*t ) * ( V \ x0 )  % THIS IS WRONG, FIX YOUR NOTES

% Most students (trusting their professors as they probably shouldn't!)
% used the wrong command above in their HW and got almost the right answer.

% The problem is this:  exp( D*t ) takes every element of D*t and raises
% it to the e power.  On the diagonal this gives what is desired. But
% off the diagonal it gives  exp(0) which is  1 (one) in every position
% where what we want is zero.

% Why does the wrong command give  almost the right answer?
% Because in the assigned problem the diagonal elements of 
% exp(D*t) were so big that the  ones off the diagonal were effectively
% almost like the desired zeros.

% A command, almost looking like the wrong command above, but that happens
% to be right, is:

x_of_t = V  *   expm( D*t ) * ( V \ x0 )  %THIS IS CORRECT

% Here the command EXPM is used.  It is the matrix exponential, something
% you can read about in EDWARDS and PENNEY section 5.5 (which we didn't
% cover). The matrix exponential is NOT the matrix each of whose
% components are raised to the e power.  In the case of a diagonal 
% matrix (like D or D*t) the matrix exponential IS the diagonal
% components raised to the e power with the zero's off the diagonal
% being left alone.

% Why did Andy write that wrong command and then convince Professor
% Terrell, TA Tamas, and hundreds of students to do something wrong?
% Because he wrote out the correct command  with DIAG above (thats
% why he put DIAG in the commands in the MATLAB matrix examples file)
% and then, on a later day under the influence of a stray cosmic
% ray, 'realized' (incorrectly) that exp(0)=0.  This led to the 
% erronious conclusion that the two DIAG commands could be eliminated.