% Programming and Numerical Methods
% MSc in Computational and Applied Physics

% Code 23: Explicit Runge-Kutta 4 ODE solver
% A. Meseguer (2007)

% On input: fut (vector field to be called f(z))
%               z = [u_1 ; u_2 ; ... u_N ; t]
%           u0 (initial condition u0(N,1))
%           k  (time step)
%           T  (integration time) 
% On output: t (array of time grid)
%            u (array of v_k's)


  function [t,u] = rk4(fut,u0,k,T)

  t = [0] ; u = [u0] ; time = 0 ;

  while time < T
   z = [u(:,end) ; t(end)] ;
   a = k*feval(fut,z) ; 
   z = [u(:,end)+a/2 ; t(end)+k/2] ;
   b = k*feval(fut,z) ; 
   z = [u(:,end)+b/2 ; t(end)+k/2] ;
   c = k*feval(fut,z) ; 
   z = [u(:,end)+c ; t(end)+k] ;
   d = k*feval(fut,z) ; 
   u = [u , u(:,end) + (a+2*b+2*c+d)/6] ;
   time = time + k ; t = [t time] ;
  end  

  
  
