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

% x, y = initial values; y must be a row vector.
% xStop = terminal value of x.
% H = increment of x at which solution is stored.
% tol = per-step error tolerance (default = 1.0e-6).
% OUTPUT:
% xSol, ySol = solution at increments H.

ifsize(y,1)>1;y=y';end %ymustberowvector
if nargin < 6; tol = 1.0e-6; end
n = length(y);
xSol = zeros(2,1); ySol = zeros(2,n);
xSol(1) = x; ySol(1,:) = y;
k=1;
whilex<xStop
k=k+1;
H = min(H,xStop - x);
y = midpoint(dEqs,x,y,x + H,tol);
x=x+H;
xSol(k) = x; ySol(k,:) = y;

end

EXAMPLE 7.10
Compute the solution of the initial value problem

y = sin y y (0) = 1

at x 4, followedbyRichardson

extrapolation (this problemwas solvedwith the second-order Runge-Kuttamethod
in Example 7.3).

= 0 . 5with the midpointformulas using n = 2 and n =

Solution With n = 2 the step length is h = 0 . 25. The midpointformulas, Eqs. (7.26)

and (7.27), yield

y 1 = y 0 + hf 0 = 1 + 0 . 25 sin 1 . 0 = 1 . 210 368

y 2 = y 0 + 2 hf 1 = 1 + 2(0 . 25)sin 1 . 210 368 = 1 . 467 87 3

1
2 ( y 1 +

. =

y 0 +

y h (0 5) hf 2 )

1
2 (1

= . + . + . .

210 368 1 467 87 3 0 25 sin 1 467 87 3)

= 1 . 463 459