function [T,Y] = rk4(f,a,b,ya,m)
% [T,Y] = rk4(f,a,b,ya,m)
% Runge-Kutta solution for y' = f(t,y) with y(a) = ya.
% f  is the function, input.
% a  is the left endpoint, input.
% b  is the right endpoint, input.
% ya is the initial condition, input.
% m  is the number of steps, input.
% T  is the vector of abscissas, output.
% Y  is the vector of ordinates, output.
h = (b - a)/m;
T = zeros(1,m+1);
Y = zeros(1,m+1);
T(1) = a;
Y(1) = ya;
for j=1:m,
  tj = T(j);
  yj = Y(j);
  k1 = h*feval(f,tj,yj);
  k2 = h*feval(f,tj+h/2,yj+k1/2);
  k3 = h*feval(f,tj+h/2,yj+k2/2);
  k4 = h*feval(f,tj+h,yj+k3);
  Y(j+1) = yj + (k1 + 2*k2 + 2*k3 + k4)/6;
  T(j+1) = a + h*j;
end
