% ADAMS-Bathforth FOUR-STEP  ALGORITHM
 %
 % 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 FOUR-STEP 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 4th-order 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 AB METHOD\n\n');
 
% STEP 6, Using AB 4-step method to compute approximate solution
 for I = 4:N 
% STEP 7

 T(I+1) = T(I)+H;

 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)) - 9.0*F(T(I-3),W(I-3));
 W(I+1) = W(I)+H*(Part1)/24.0;
 fprintf(OUP, '%2d   %5.3f %11.7f\n', I+1, T(I+1), W(I+1));
 
 end;