echo off;
% NUMERICAL METHODS: MATLAB Programs, (c) John H. Mathews 1994
% To accompany the text:
% NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992
% Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A.
% This free software is complements of the author.

% Algorithm 9.5 (Runge-Kutta-Fehlberg Method RKF45).
% Section 9.5, Runge-Kutta Methods, Page 461
echo on; clc; format long; hold off; clear

% This program implements Runge-Kutta-Fehlberg method

% for solving the initial value problem

% y' = f(t,y) with y(a) = ya

% Define and store the function f(t,y) in the M-file f.m
% function z = f(t,y)
% z = 1 + y^2;

delete f.m
diary f.m; disp('function z = f(t,y)');...
disp('z = 1 + y^2;');...
diary off;

f(0,0); % Remark. f.m and rkf45.m are used for Algorithm 9.5

pause % Press any key to continue.

clc;
% Place the endpoints of [a,b] in a and b.

% Place the initial value y(a) in ya.

% Place the tolerance in toler.

% Place the tentative number of subintervals in m.

a = 0;
b = 1.4;
ya = 0;
toler = 1e-5;
m = 14;

[T,Y] = rkf45('f',a,b,ya,m,toler);
points = [T;Y];

pause % Press any key to continue.

clc;, clg;
c = 0;
d = 6;
axis([a b c d]);...
plot(T,Y,'g');...
hold on;...
plot([a b],[0 0],'b',[0 0],[c d],'b');...
if m<=30,
plot(T,Y,'or');
end;...
xlabel('t');...
ylabel('y');...
title('Runge-Kutta-Fehlberg solution to y` = f(t,y)');...
grid;...
hold off;...
shg; pause % Press any key to continue.

Mx1 = 'Runge-Kutta-Fehlberg solution to y` = f(t,y).';
Mx2 = ' t(k) y(k)';
clc,echo off,diary output,...
disp(''),disp(Mx1),...
disp(''),disp(Mx2),disp(points'),diary off,echo on