echo off;
% NUMERICAL METHODS by John H. Mathews
% Algorithm 9.4 (Runge-Kutta Method of Order 4).
% Section  9.5,  Runge-Kutta Methods, Page 460
% PRENTICE HALL, INC. © 1992
echo on; clc; format long; hold off; clear

% This program implements the Runge-Kutta method

% for solving the initial value problem

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

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

delete f.m
diary f.m; disp('function z = f(t,y)');...
           disp('z = (t-y)/2;');...
diary off;

f(0,0); % Remark. f.m and rk4.m are used for Algorithm 9.4
pause   % Press any key to continue.

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

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

%  Place the number of subintervals in  m.

a  = 0;
b  = 3;
ya = 1;
m  = 6;

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

%pause  % Press any key to see the list of solution points.


kuvano=1;
figure(kuvano);
clc;, clg;
c = 0;
d = 1.75;
axis([a b c d]);...
plot(T,Y,'-g');...
hold on;...
plot([a b],[0 0],'b',[0 0],[c d],'b');...

if m<=32,
  plot(T,Y,'or');
end;...
xlabel('t');...
ylabel('y');...
title('Runge-Kutta solution to y` = f(t,y)');...
grid;...
axis;...
hold off;...
shg; pause  % Press any key to continue.

Mx1 = 'Runge-Kutta 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

pause % press any ..., sitten voit komentaa: rk4more
