% x' = A x


function rk2_driver()


   A = [ -21  19 -20; 19 -21 20; 40 -40 -40];
   dt = 0.02;
   finalTime = 1;
   T = 0:dt:finalTime;
   v = zeros(3,length(T));
   v(:,1) = [1; 0; -1];    % initial condition
   ex = zeros(3,length(T));
   ex(:,1) = [1; 0; -1]; 

   t = 0; k = 1;
   while t<finalTime

     v(:,k+1) = rk2(v(:,k),t,dt,@F);
     ex(:,k+1) = exact(t+dt);
  
     t = t + dt;
     k = k + 1;
      
   end
   
   
   
   plot(T,v(1,:),'b',T,ex(1,:),'r--')


  function fp = F(u,t)
      fp = A*u;
  end

  function f = exact(t)

     f(1) =  0.5*exp(-2*t) + 0.5*exp(-40*t)*(cos(40*t) + sin(40*t));
     f(2) =  0.5*exp(-2*t) - 0.5*exp(-40*t)*(cos(40*t) + sin(40*t));
     f(3) = -exp(-40*t)*( cos(40*t)-sin(40*t) );
     f = f(:);

  end

end