
function normalmodes ()

%Parameters
m1=1; m2=1;
k01 = 1; k12=1;
M = [ m1 0; 0 m2];
K = [ (k01+k12)  -k12; -k12 k12];

%ICs:
%x0= [1 1.618]'; % initial position
x0= [1 -.618]'; % initial position
v0=[ 0 0]'; % initial velocity
zzero= [x0;v0];

n= 1000;

tspan = linspace(0, 20, n);

%Method I: ODE soln
[t zmatrix] = ode45(@rhs, tspan, zzero, [], M,K );
  x1 = zmatrix(:,1);
  x2 = zmatrix(:,2);

  clf
  figure(1)

  plot (t,x1,'r.',t,x2,'b');  title('free response')
  
  
%Math method I: using M^(-1).
zmatrix = math(tspan, x0,v0, M,K);
    x1 = zmatrix(:,1);
    x2 = zmatrix(:,2);

  figure(2)

  plot (t,x1,'g.',t,x2,'k');  title('free response')
  
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function zdot = rhs(t,z,M,K)
x = z(1:2);
v = z(3:4);

xdot = v;

vdot = -M^-1 * K*x;


zdot = [xdot; vdot];

end
%%%%%%%%%%

%%%%%%%%%%
function zmatrix = math(tspan, x0,v0, M,K)

%Find e-vectors and e-values of M^(-1)
[ P omsq] = eig(M^(-1)*K);

%For a simple example, assume v0=0
  A = P\x0;         %Solve for ICs
  n= length(tspan);
  zmatrix= zeros(n,4);
  
  %Evaluate positions at each time
  for i=1:n
  
   cosomt= cos((diag(omsq)).^(1/2)*tspan(i)) ; 
      
      x= P*(A.*cosomt); %Positions at time t
  
  zmatrix(i,[1 2]) = x'; %First two entries are positions
  end

end
%%%%%%%%%%%%%