% ADAMS-FOURTH ORDER PREDICTOR-CORRECTOR ALGORITHM 5.4
 %
 % To approximate the solution of the initial value problem
 %        y' = f(t,y), a <= t <= b, y(a) = alpha,
 % at N+1 equally spaced points in the interval [a,b].
 %
 % INPUT:   endpoints a,b; initial condition alpha; integer N.
 %
 % OUTPUT:  approximation w to y at the (N+1) values of t.
 syms('F', 'OK', 'A', 'B', 'ALPHA', 'N', 'FLAG', 'NAME', 'OUP');
 syms('H', 'T', 'W', 'I', 'K1', 'K2', 'K3', 'K4', 'T0', 'W0', 'J');
 syms('t','y', 's','Part1','Part2');
 
 s = 'y-t^2+1';
 F = inline(s,'t','y');
 
 A = 0.0; % starting time
 B = 2;  %ending time
 
 ALPHA = 0.5; %Initial condition
 

 N = 10;  % number of subintervals
 T = zeros(1,N+1);
 W = zeros(1,N+1);
 
 OUP = 1;

 fprintf(OUP, 'ADAMS-BASHFORTH FOURTH ORDER PREDICTOR CORRECTOR METHOD\n\n');
 fprintf(OUP, 'step    t           w\n');
% STEP 1
 H = (B-A)/N;
 T(1) = A;
 W(1) = ALPHA;
 fprintf(OUP, '%2d   %5.3f %11.7f\n', 1, T(1), W(1));
 
% STEP 2
 for I = 1:3 
% STEP 3 AND 4
% compute starting values using Runge-Kutta method
 T(I+1) = T(I)+H;
 K1 = H*F(T(I), W(I));
 K2 = H*F(T(I)+0.5*H, W(I)+0.5*K1);
 K3 = H*F(T(I)+0.5*H, W(I)+0.5*K2);
 K4 = H*F(T(I+1), W(I)+K3);
 W(I+1) = W(I)+(K1+2.0*(K2+K3)+K4)/6.0;
% STEP 5
 fprintf(OUP, '%2d   %5.3f %11.7f\n', I+1, T(I+1), W(I+1));
 end;
 
 fprintf(OUP, 'Start with PREDICTOR CORRECTOR METHOD\n\n');
 
% STEP 6
 for I = 4:N 
% STEP 7
% W0 will be used in place of t, w resp.
 
 T(I+1) = T(I)+H;
 
% predict W(I+1)
 Part1 = 55.0*F(T(I),W(I))-59.0*F(T(I-1),W(I-1))+37.0*F(T(I-2),W(I-2));
 Part2 = -9.0*F(T(I-3),W(I-3));
 W0 = W(I)+H*(Part1+Part2)/24.0;
 
% correct W(I+1)
 Part1 = 9.0*F(T(I+1),W0)+19.0*F(T(I),W(I))-5.0*F(T(I-1),W(I-1))+F(T(I-2),W(I-2));
 W(I+1) = W(I)+H*(Part1)/24.0;
% STEP 8
 fprintf(OUP, '%2d   %5.3f %11.7f\n', I+1, T(I+1), W(I+1));
 
 end;