clear;									% again, clear all variables
%
% we want to do a similar test, so make a list of timestep sizes to look at
%
hlist=[0.05 0.025 0.0125];
%
% declare the rate function f as an inline funciton, similar to 5.1
% we also know the exact solution (by hand or by maple) so type that in too
%
f=inline('2*y-3*t','t','y');
yexact=inline('3*t/2 + 3/4 + (1/4)*exp(2*t)','t');
%
% Now start the main loop.  Loop through
% 3 timestep sizes listed in hlist
%
for ih=1:3
	h=hlist(ih);			% let h be the current h being tested
	disp(h);				% display it
	t=0:h:1;				% make a vector of t values
	N=length(t);			% again it is of length N
	% since we're using AM4 we will be starting the
	% timestepping from t(5), meaning we should write
	% down yex for the first 4 times since it will not
	% be covered in the subsequent loop
	yex(1)=yexact(t(1));	% Exact
	yex(2)=yexact(t(2));	% Exact
	yex(3)=yexact(t(3));	% Exact
	yex(4)=yexact(t(4));	% Exact
	% 
	% now for the approximate solution, y
	%
	% since we're using AM4 and AB4, we cannot start
	% from t(2) (the second timestep) since these methods
	% use information from the previous 4 timestep.  This
	% should be more accuarate, right?  We shall see...
	%
	% So use RK4 to get the first 4 approximate values
	% i.e.  find y at times t(2) t(3) and t(4);
	%
	% initial condition
	y(1)=1;
	% No RK4
	k1=f(t(1),y(1)); k2=f(t(1)+0.5*h,y(1)+0.5*h*k1); k3=f(t(1)+0.5*h,y(1)+0.5*h*k2); k4=f(t(1)+h,y(1)+h*k3);
	y(2)=y(1)+(h/6)*(k1+2*k2+2*k3+k4);
	k1=f(t(2),y(2)); k2=f(t(2)+0.5*h,y(2)+0.5*h*k1); k3=f(t(2)+0.5*h,y(2)+0.5*h*k2); k4=f(t(2)+h,y(2)+h*k3);
	y(3)=y(2)+(h/6)*(k1+2*k2+2*k3+k4);
	k1=f(t(3),y(3)); k2=f(t(3)+0.5*h,y(3)+0.5*h*k1); k3=f(t(3)+0.5*h,y(3)+0.5*h*k2); k4=f(t(3)+h,y(3)+h*k3);
	y(4)=y(3)+(h/6)*(k1+2*k2+2*k3+k4);
	% so y is known through 4 timesteps.  We can now begin the AM4/AB4
	% timestepping
	for n=4:(N-1)
	% Euler Predictor (comment out if not used)
	%	y(n+1) = y(n) + h*f(t(n),y(n));
	%
	%
	% AB4 Predictor (comment out if not used)
	% notice this follows the book exactly and
	% we will store the predicted value at t(n+1)
	% as y(n+1)
		y(n+1) = y(n) + (h/24)*(55*f(t(n),y(n))...
								-59*f(t(n-1),y(n-1))...
								+37*f(t(n-2),y(n-2))...
								 -9*f(t(n-3),y(n-3)));
		%correct=0;		% no correction steps
		correct=1;		% 1 correction step
		%
		% loop through the correction steps using AM4
		% this also exactly follows the formula in the book
		% note: we are always overwriting y(n+1) so the predition
		% y(n+1) is no longer available.
		%
		for ci=1:correct
			y(n+1) = y(n) + (h/24)*(9*f(t(n+1),y(n+1))...
									+19*f(t(n),y(n))...
									-5*f(t(n-1),y(n-1))...
									+f(t(n-2),y(n-2)));
		end
		yex(n+1) = yexact(t(n+1));
	end
	% list the error for each h value at t=1.0 (this is at the 
	% end of y and yex
	%
	% notice no ;  suppressing the semicolon causing output to the matlab
	% window
	y(end)-yex(end)
	%
	% now plot
	%
	figure;					% open a new figure
	plot(t,y-yex);  % plot t versus y-yexact=global error
					% does the error go down when h is decreased?
end