% Program 6.7 Multistep method plotting program
% Inputs: [inter(1),inter(2)] time interval,
%  ic = [y0]  initial condition,
%  h=stepsize, s=number of (multi)steps, e.g. 2 for 2-step method
% Output: solution y
% Calls a multistep method such as ab2step.m
% Example usage: y=exmultistep([0 1],1,0.05,2);
function [t,y]=exmultistep(inter,ic,h,s)
n=round((inter(2)-inter(1))/h);
% Start-up phase
y(1,:)=ic;t(1)=inter(1);
for i=1:s-1               % start-up phase, using one-step method
  t(i+1)=t(i)+h;
  y(i+1,:)=trapstep(t(i),y(i,:),h);
  f(i,:)=ydot(t(i),y(i,:));
end
for i=s:n                 % multistep method loop
  t(i+1)=t(i)+h;
  f(i,:)=ydot(t(i),y(i,:));
  y(i+1,:)=unstable2step(t(i),i,y,f,h);
end

function y=trapstep(t,x,h)
%one step of the Trapezoid Method from section 6.2
z1=ydot(t,x); 
g=x+h*z1; 
z2=ydot(t+h,g); 
y=x+h*(z1+z2)/2;

function z=ab2step(t,i,y,f,h)
%one step of the Adams-Bashforth 2-step method
z=y(i,:)+h*(3*f(i,:)/2-f(i-1,:)/2);

function z=unstable2step(t,i,y,f,h)
%one step of an unstable 2-step method
z=-y(i,:)+2*y(i-1,:)+h*(5*f(i,:)/2+f(i-1,:)/2);

function z=weaklystable2step(t,i,y,f,h)
%one step of a  weakly-stable 2-step method
z=y(i-1,:)+h*2*f(i,:);

function z = ydot(t,y)  % IVP from section 6.1
z = t*y+t^3;
