% Here is an example for RK4 for solving the first order system of ODEs
% in B&F 8th ed, 5.9.1 a
% which is given by
% 
%   y_1' = 3*y_1 + 2*y_2 - (2*t^2+1)*exp(2*t)
%   y_2' = 4*y_1 +   y_2 + (t^2+2*t-4)*exp(2*t)
% 
% with initial condition
% 
%   y_1(0) = 1
%   y_2(0) = 1
% 
% for 0<=t<=1 and with h=0.2. The true solution is:
%
%   y_1(t) = exp(5*t)/3 -   exp(-t)/3 + exp(2*t)
%   y_2(t) = exp(5*t)/3 + 2*exp(-t)/3 + t^2*exp(2*t)
%
% Sample output from this program (identical to solution in book)
%
%   >> rk4_example
%    t_i       y1 calc       y1 true       y2 calc       y2 true
%    0.0    1.00000000    1.00000000    1.00000000    1.00000000
%    0.2    2.12036583    2.12500839    1.50699185    1.51158743
%    0.4    4.44122776    4.46511961    3.24224021    3.26598528
%    0.6    9.73913329    9.83235869    8.16341700    8.25629549
%    0.8   22.67655977   23.00263945   21.34352778   21.66887674
%    1.0   55.66118088   56.73748265   56.03050296   57.10536209


% setup function f such that y' = f(t,y) (where y is a vector)
f = @(t,y) [ 3*y(1) + 2*y(2) - (2*t^2+1)*exp(2*t)
             4*y(1) +   y(2) + (t^2+2*t-4)*exp(2*t) ];

% setup initial condition y(a) = y0
y0 = [1;1];

% use 5 points to have h=0.2
n=5; 

% call Runge-Kutte 4 for systems code
y = rk4(0,1,n,y0,f);

% true solution
y_1 = @(t) exp(5*t)/3 -   exp(-t)/3 + exp(2*t);
y_2 = @(t) exp(5*t)/3 + 2*exp(-t)/3 + t^2*exp(2*t);

% print result as a table
fprintf('%4s %13s %13s %13s %13s\n','t_i','y1 calc','y1 true','y2 calc','y2 true');
for i=1:n+1,
 ti = (i-1)*h;
 fprintf('%4.1f %13.8f %13.8f %13.8f %13.8f\n',ti,y(1,i),y_1(ti),y(2,i),y_2(ti));
end;


