% the program takes the initial data y_0 at the initial time x_0 and 
% computes up to time x_final using N timesteps of size dx.  It creates
% two vectors, x and y, which you can then plot against each other
% plot(x,y)
%
% To call the routine, use the command:
%
% [x,y] = runge_kutta_4(y_0,x_0,x_final,N,@f)
%
function [x,y] = runge_kutta_4(y_0,x_0,x_final,N,f)

% step-size:
dx = (x_final-x_0)/N;

% for dy/dx = f(x,y(x))
%
% first define 
%         F1 = f(x,y)
%         F2 = f(x+dx/2,y+dx/2*F1)
%         F3 = f(x+dx/2,y+dx/2*F2)
%         F4 = f(x+dx,y+dx*F3)
% y(x+dx) ~ y(x) + dx/6*(F1 + 2*F2 + 2*F3 + F4)

% pux initial data into vector
x(1) = x_0;
y(:,1) = y_0;
for ii=1:N
  F1 = feval(f,x(ii),y(:,ii));
  F2 = feval(f,x(ii)+dx/2,y(:,ii)+dx/2*F1);
  F3 = feval(f,x(ii)+dx/2,y(:,ii)+dx/2*F2);
  F4 = feval(f,x(ii)+dx,y(:,ii)+dx*F3);
  y(:,ii+1) = y(:,ii) + dx/6*(F1 +  2*F2 + 2*F3 + F4);
  x(ii+1) = x_0 + ii*dx;
end