% This code tests whether the order  n  of the starting method 
%  for the multistep method of order  m  can satisfy
%          n = m-1
%  and still not degrade the order of the multistep method.

% I'll do it for the 3rd-order Admas--Bashforth method 
%        y_{i+1} - y_i = (h/12)*(23*f_i - 16*f_{i-1} + 5 f_{i-2})
%  using the modified Euler, midpoint, and simple Euler methods as starting methods
%  to obtain  y_1  and  y_2  and hence  f_1 and f_2.


global a           % This makes parameter  a  available to all functions.

% Set the parameters of the problem:
a=1;

x0=0; 
xf=1;                          % beginning and ending values of x
h=0.025;                         % step size in x
N=round((xf-x0)/h);     % total # of integration steps
xf_starting=2*h;               % ending value of x for the starting method

y0=1;                          % initial value of y

xspan_starting=[x0 xf_starting];    % span for the starting method
xspan = [x0 xf];                    % whole span


% Solve the equation using the starting method:
[x_sm, y_sm]=ME_1eq(@fun4_hw3_starting_method,y0,xspan_starting,h);  % mod Euler
% [x_sm, y_sm]=midpoint(@fun4_hw3_starting_method,y0,xspan_starting,h);   % midpoint

% [x_sm, y_sm]=hw1_EX_Euler_callsfun(@fun4_hw3_starting_method,xspan_starting,h,y0);  % simple Euler

                               % This invokes the function which runs one of one-step methods
                               %  to integrate the function supplied in the function file  
                               %  fun4_hw3_starting_method.m
                               
[x, Y]=hw3_AdamsBashforth3(@fun4_hw3_starting_method,xspan,h,y0,y_sm(2),y_sm(3));

% y_exact = ((x+2)/2).^2;   % exact solution of y'=sqrt(y),  y(0)=1
% y_exact = tan((a/2)*x.^2+atan(y0));  % exact solution of y'=ax(1+y^2),  y(0)=1
% y_exact = atan(a*x.^2+tan(y0));  % exact solution of y'=ax(1+cos(2y)),  y(0)=1
% y_exact = tan((a/4)*x.^4+atan(y0));  % exact solution of y'=ax^3(1+y^2),  y(0)=1
y_exact = -4*exp(2*x)+5*exp(3*x);  % exact solution of y'=2y+5*exp(3x),  y(0)=1



figure(13);
plot(x,Y,'mo','markersize',5)
hold on
plot(x,y_exact,'k')
hold off

GE=y_exact(end)-Y(end)