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

k=1;
whilex<xStop
h = min(h,xStop - x);
d = feval(deriv,x,y); % Derivatives of [y]

hh=1;
forj=1:4 %BuildTaylorseries

hh = hh*h/j; % hh = hˆj/j!

y=y+d(j,:)*hh;

end
x=x+h;k=k+1;
xSol(k) = x; ySol(k,:) = y; % Store current soln.

end

printSol

Thisfunctionprints the results xSol and ySol in tabular form. The amountofdata
iscontrolledbythe printoutfrequency freq .For example, if freq = 5 ,every fifth
integration stepwouldbe displayed. If freq = 0 ,only the initial and final values will
be shown.

function printSol(xSol,ySol,freq)
% Prints xSol and ySoln arrays in tabular format.
% USAGE: printSol(xSol,ySol,freq)
% freq = printout frequency (prints every freq-th
%

line of xSol and ySol).

[m,n] = size(ySol);
if freq == 0;freq = m; end
head = '

x';

fori=1:n
head = strcat(head,' y',num2str(i));

end
fprintf(head); fprintf('\n')
fori=1:freq:m
fprintf('%14.4e',xSol(i),ySol(i,:)); fprintf('\n')

end
if i ˜= m; fprintf('%14.4e',xSol(m),ySol(m,:)); end

EXAMPLE 7.1
Given that

y + x 2 4 y = y (0) = 1