Initial Value Problems - Numerical Methods in Engineering with MATLAB - page 287

% OUTPUT:
%y=y(xStop).
ifsize(y,1)>1;y=y';end %ymustberowvector
if nargin <5; tol = 1.0e-6; end
kMax = 51;
n = length(y);
r = zeros(kMax,n); % Storage for Richardson extrapolation.
% Start with two integration steps.
nSteps = 2;
r(1,1:n) = mid(dEqs,x,y,xStop,nSteps);
rOld = r(1,1:n);
fork=2:kMax
% Double the number of steps & refine results by
% Richardson extrapolation.
nSteps = 2*k;
r(k,1:n) = mid(dEqs,x,y,xStop,nSteps);
r = richardson(r,k,n);
% Check for convergence.
dr = r(1,1:n) - rOld;
e = sqrt(dot(dr,dr)/n);
ife<tol;y=r(1,1:n);return;end
rOld = r(1,1:n);
end
error('Midpoint method did not converge')
functionr=richardson(r,k,n)
% Richardson extrapolation.
forj=k-1:-1:1
c =(k/(k-1))ˆ(2*(k-j));
r(j,1:n) =(c*r(j+1,1:n) - r(j,1:n))/(c - 1.0);
end
return
functiony=mid(dEqs,x,y,xStop,nSteps)
% Midpoint formulas.
h = (xStop - x)/nSteps;
y0=y;
y1 = y0 + h*feval(dEqs,x,y0);
fori=1:nSteps-1