function [y,t] = rk4(f,t0,y0,T,n)
% by Shelvean Kapita
% Runge-Kutta method of order 3
% y' = f(t,y)
% t0 = initial time
% y0 = initial value
% T  = final time
% n = number of time steps
h = (T-t0)/n; % time step
i = 1;
t = linspace(t0,T,n+1);
y = zeros(size(t));
y(1)=y0;
while (i<=n)
    k1 = f(t(i),y(i));
    k2 = f(t(i)+h/2, y(i)+0.5*h*k1);
    k3 = f(t(i)+h/2, y(i)+0.5*h*k2);
    k4 = f(t(i)+h, y(i)+h*k3);
    y(i+1)=y(i)+h*(k1+2*k2+2*k3+k4)/6;
    i = i+1;
end