%
%	fourth_rk.m:  runge-kutta timestep
%
%
function [rkZn,rktn] = ode_rk(rkZ,rkt,rkDt,par)

k1 = rkDt*deriv(rkZ       ,rkt         ,par);
k2 = rkDt*deriv(rkZ+0.5*k1,rkt+0.5*rkDt,par);
k3 = rkDt*deriv(rkZ+0.5*k2,rkt+0.5*rkDt,par);
k4 = rkDt*deriv(rkZ+    k3,rkt+    rkDt,par);

rkZn = rkZ + (k1 + 2*k2 + 2*k3 + k4)/6;
rktn = rkt + rkDt;

%- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
%
%  4th-order front with newton
%
function rkP = deriv(rkZ,rkt,par)

rkP = zeros(size(rkZ));

rkP(1) = rkZ(2);
rkP(2) = rkZ(3); 
rkP(3) = -0.5*(rkZ(1)^2 - par(1));

rkP(4) = rkZ(5);
rkP(5) = rkZ(6); 
rkP(6) = -rkZ(1)*rkZ(4);