%      PROGRAM IMPLICIT_EULER.M
%
  clear all

% 
% Define variables
%

  m = 4;
  k = 16;
  c = 2;

  dt = input('dt = ');    
  T = 0:dt:20;
  N = length(T);

  y0 = 2;
  y1 = 30;

  nu = sqrt(4*k*m - c^2)/(2*m);
  a1 = y0;
  a2 = (y1 + c*y0/(2*m))/nu;
  
%
% Create system matrices and vectors
%

  A = [0 1; -k/m -c/m];
  AA = inv(eye(2,2) - dt*A);
  z0 = [y0;y1];
  
%
% Iterate to approximate the solution using the implicit Euler
% algorithm
%

  z(:,1) = z0;
  y(1) = y0;
  for j=2:N
     t = T(j);
     y(j) = exp(-c*t/(2*m))*(a1*cos(nu*t) + a2*sin(nu*t));   % True solution
     z(:,j) = AA*z(:,j-1);                                   % Approximate solution
  end
  
%
% Plot the results
%

  figure(1)
  Displacement = z(1,:);
  plot(T,Displacement,'-g',T,y,'--r',T,0*Displacement,'linewidth',2);
  h = gca;
  set(h,'FontSize',[18]);
  xlabel('Time (s)')
  ylabel('Displacement (m)')
  legend('Implicit Euler','True',4)