% Programming and Numerical Methods
% MSc in Computational and Applied Physics

% Code 24: Comparing performances of different ODE solvers
% A. Meseguer (2007)

  clear all
  close all

  errexpeu = [] ; errmodeu = [] ; errheun = [] ; errrk4 = [] ; 
  kvec = [] ; T = 20 ;

  for k = [.5 .1 .05 .01 .005 .001]
   
   [t,uexpleuler] = expleuler(@fut,[1 ; 0],k,T);
   [t,umodeuler] = modeuler(@fut,[1 ; 0],k,T);
   [t,uheun] = heun(@fut,[1 ; 0],k,T);
   [t,urk4] = rk4(@fut,[1 ; 0],k,T);


   kvec = [kvec k] ; 
   errexpeu = [errexpeu abs(uexpleuler(1,end)-cos(T))] ;
   errmodeu = [errmodeu abs(umodeuler(1,end)-cos(T))] ;
   errheun = [errheun abs(uheun(1,end)-cos(T))] ;
   errrk4 = [errrk4 abs(urk4(1,end)-cos(T))] ;

  end

  loglog(kvec,errexpeu); hold on ;  loglog(kvec,errexpeu,'ob') ; grid on
  loglog(kvec,errmodeu,'-r');  loglog(kvec,errmodeu,'or')
  loglog(kvec,errheun,'-m') ;  loglog(kvec,errheun,'om') 
  loglog(kvec,errrk4,'-k')  ;  loglog(kvec,errrk4,'ok')

