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

forp=2:5000
% Compute K's from Eq. (7.18)
K = zeros(6,n);
K(1,:) = h*feval(dEqs,x,y);
fori=2:6
BK = zeros(1,n);
forj=1:i-1
BK = BK + B(i,j)*K(j,:);
end
K(i,:)=h*feval(dEqs,x+A(i)*h,y+BK);
end
% Compute change in y and per-step error from
% Eqs.(7.19) & (7.20)
dy=zeros(1,n);E=zeros(1,n);
fori=1:6
dy = dy + C(i)*K(i,:);
E=E+(C(i)-D(i))*K(i,:);
end
e = sqrt(sum(E.*E)/n);
% If error within tolerance, accept results and
% check for termination
ife<=eTol
y=y+dy;x=x+h;
k=k+1;
xSol(k) = x; ySol(k,:) = y;
if stopper == 1;
break
end
end
% Size of next integration step from Eq. (7.24)
if e˜= 0; hNext = 0.9*h*(eTol/e)ˆ0.2;
else; hNext=h;
end
% Check if next step is the last one (works
% with positive and negative h)
if(h>0)==(x+hNext>=xStop)
hNext = xStop - x; stopper = 1;
end
h=hNext;
end