function ynew = ruka4(deriv,y,h,t)
%RUKA4  4th order Runge-Kutta solver for y' = f(t,y) 

%   YNEW = RUKA4(DERIV,Y,H,T) advances the solution vector Y over an
%   interval H and returns the incremented variables YNEW.
%
%   User supplies the function DERIV(T,Y) which returns the 
%   derivatives DERIV at time T.
%  
%   See also DERIVS, SIGMNETWORK.
%

%   Copyright 2003
%   Centre for Nonlinear Dynamics in Physiology and Medicine

hh = h*0.5;
h6 = h/6;
th = t+hh;

dydt = feval(deriv,t,y);
yt = y+hh*dydt;

dyt = feval(deriv,th,yt);
yt = y+hh*dyt;

dym = feval(deriv,th,yt);
yt = y+h*dym;
dym = dyt+dym;

dyt = feval(deriv,t+h,yt);
ynew = y+h6*(dydt+dyt+2*dym);

return



