% Igor Yanovsky
% Teaching Assistant
% Math 151B
% 04/22/07
%
%#####################################################################
%
% This program solves the ordinary differential equation
% 
%  dy/dt = f(t,y)
%  y(t_0) = y_0
% 
% using ADAMS FOURTH-ORDER PREDICTOR-CORRECTOR method.
% The function defining the ODE is defined in a separate Matlab function.
%
%#####################################################################

clear all

format long

a =  0.0;
b =  1.0;
y0 = 1.0;

h = 0.1;

N = (b - a)/h;

t = a;
w(1) = y0;
yexact(1) = y0;
error(1) = 0;

for(i=1:3)
    % advances the ODE using RK-4:
    k1 = h*f( t        , w(i)          );
    k2 = h*f( t+(h/2.0), w(i)+(k1/2.0) );
    k3 = h*f( t+(h/2.0), w(i)+(k2/2.0) );
    k4 = h*f( t+h      , w(i)+k3       );
    
    w(i+1) = w(i) + (1.0/6.0)*( k1 + 2.0*k2 + 2.0*k3 + k4 );
    
    t = t + h;
 
    yexact(i+1) = findExact( t );
    error(i+1) = abs( w(i+1)-yexact(i+1) );
end

% Adams Fourth-Order Predictor-Corrector method:

for(i=4:N)
    % Adams-Bashforth predictor step:
    w_tmp = w(i) + (h/24) * ( 55*f(t,w(i)) - 59*f(t-h,w(i-1)) + 37*f(t-2*h,w(i-2)) - 9*f(t-3*h,w(i-3)) );
    
    % Adams-Moulton corrector step:
    w(i+1) = w(i) + (h/24) * ( 9*f(t+h,w_tmp) + 19*f(t,w(i)) - 5*f(t-h,w(i-1)) + f(t-2*h,w(i-2)) );
    
    t = t + h;
    
    yexact(i+1) = findExact( t );
    error(i+1) = abs( w(i+1)-yexact(i+1) );
end


fprintf(1,'N = %3d',N);  fprintf(1,', h = %4.7e',h);  fprintf(1,', t = %4.2f',t);  fprintf(1,',\n');

w
yexact
error


fprintf(1,'w = %5.10e', w(N+1)    );  fprintf(1,',\n');
fprintf(1,'y = %5.10e', yexact(N+1));  fprintf(1,',\n');
fprintf(1,'error = %5.10e', error(N+1) );  fprintf(1, '\n');


clear all