%      PROGRAM BACK_EULER.M
%
  clear all
  close all

% 
% Define variables
%

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

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

  x0 = 2;
  x1 = 30;

  nu = sqrt(4*k*m - c^2)/(2*m);
  a1 = x0;
  a2 = (x1 + c*x0/(2*m))/nu;
  
%
% Create system matrices
%

  A = [0 1; -k/m -c/m];
  AA = inv(eye(2,2) - dt*A);
  X0 = [x0;x1];
  
%
% Iterate to approximate the solution using the backward Euler
% algorithm
%

  X(:,1) = X0;
  x(1) = x0;
  for j=2:N
     t = T(j);
     x(j) = exp(-c*t/(2*m))*(a1*cos(nu*t) + a2*sin(nu*t));
     X(:,j) = AA*X(:,j-1);
  end
  
%
% Plot the results
%

  figure(1)
  Displacement = X(1,:);
  plot(T,Displacement,'-g',T,x,'--r',T,0*Displacement,'linewidth',2);
  h = gca;
  get(h);
  set(h,'FontSize',[18]);
  xlabel('Time (s)')
  ylabel('Displacement (m)')
  legend('Backward Euler','True',4)
  
%  figure(2)
%  Velocity = X(2,:);
%  plot(T,Velocity,T,0*Velocity);
%  xlabel('Time (s)')
%  ylabel('Velocity (m/s)')