%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Forward Euler method
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for MTH 452-552
%%% this is a simple example how to use error estimates
%%% for simple explicit one-step methods for ODEs 
%%% using Richardson's extrapolation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%% M. Peszynska,  2005/1
%% Copyright Department of Mathematics, Oregon State University
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%% INSTRUCTIONS: calling sequence
% error_estimate_demo (flag,h)
% where flag is 1, 2, 4 (sample methods of order 1, 2, 4

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
function error_estimate_demo (flag,h)
if nargin < 2 flag = 1; h= .1; end;
%% set up the function
f = inline('3.*exp(-3*t).*y','t','y');
phi = inline('exp(1-exp(-3*t))','t');

%f = inline('-2.*t.*y.^2','t','y');
%phi = inline('1./(1+t.^2)','t');

%f   = inline('-y','t','y');
%phi = inline('exp(-t)','t');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% compute the solution and plot
%h = .05; flag = 2;

%% the numerical solution will be in y1 and y2 (double step)
[t2,y2,p] = explicit_method (f,1,0,1,2*h,flag);
n = size(t2,2);

[te,ye,p] = explicit_method (f,1,0,1,h,  flag);
t1 = zeros(n,1); y1 = zeros(n,1);

for i=1:n 
   t1(i)=te(2*i-1);y1(i)=ye(2*i-1);
end;

plot(t1,phi(t1),'k',t1,y1,'r',t2,y2,'g');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% error and estimate

%p = 1; %% Euler explicit 
%p = 2; %% any second order: modified Euler, trapezoid etc. 
%p = 4; %% Runge-Kutta: classical 

%% show error and its estimate 
for i=1:n
   fprintf('ex_error=%g err_estimate (Richardson)=%g \n',...
       y1(i) - phi(t1(i)), ...
       (y1(i) - y2(i))./(1 - 2^p) ...
);
end;

%%%%%%%%%%%%%%%%%%%%%%%
function [t,sol,p] = explicit_method (f,y0,T0,T,h,flag)

switch flag 
  case 1, [t,sol] = euler_explicit (f,y0,T0,T,h); p=1;
  case 2, [t,sol] = euler_modified (f,y0,T0,T,h); p=2; 
  case 4, [t,sol] = rk_classic (f,y0,T0,T,h);     p=4;
end;


%%%%%%%%%%%%%%%%%%%%%%%
function [t,sol] = euler_explicit (f,y0,T0,T,h)
%
%
%
%% set up the time interval and its discretization
n = floor(T-T0) / h;

%% time steps are stored in t
t = T0:h:T; 

%% the numerical solution will be in yh
sol = zeros(n+1,1);
sol (1) = y0;

%% Euler method 
for i=1:n 
   sol(i+1) = sol(i) + h * feval (f, t(i), sol(i) ); 
end;

%%%%%%%%%%%%%%%%%%%%%%%
function [t,sol] = euler_modified (f,y0,T0,T,h)
%
%
%
%% set up the time interval and its discretization
n = floor(T-T0) / h;

%% time steps are stored in t
t = T0:h:T; 

%% the numerical solution will be in yh
sol = zeros(n+1,1);
sol (1) = y0;

%% Euler modified method 
for i=1:n 
   soli = feval (f, t(i), sol(i) ); 
   sol(i+1) = sol(i) + h * feval (f, t(i)+h/2, sol(i)+h/2*soli ); 
end;


%%%%%%%%%%%%%%%%%%%%%%%
function [t,sol] = rk_classic (f,y0,T0,T,h)
%
%
%
%% set up the time interval and its discretization
n = floor(T-T0) / h;

%% time steps are stored in t
t = T0:h:T; 

%% the numerical solution will be in yh
sol = zeros(n+1,1);
sol (1) = y0;

%% classic Runge-Kutta
for i=1:n 
   soli = sol(i);
   sol1 = feval (f, t(i),       soli);
   sol2 = feval (f, t(i) + h/2, soli + h/2 * sol1); 
   sol3 = feval (f, t(i) + h/2, soli + h/2 * sol2); 
   sol4 = feval (f, t(i) + h,   soli + h   * sol3); 

   sol(i+1) = soli + h/6 * (sol1 + 2 * sol2 + 2 * sol3 + sol4);

end;