%
%	ode_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;

%- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
%
%  regularized nonlinear wave
%
function rkP = deriv(rkZ,rkt,par)

global ik

N = size(rkZ,1);
c = par(1);
b = par(2);
g = par(3);
d = par(4);

%  reconstruct ug

fact = exp(rkt*(c*ik-g*ik.^2+d*ik.^3));
vg = rkZ./fact;
ug = real(ifft(fftshift(vg))*N);

rkP = (-0.5*b) * fact .* ik .* fftshift(fft(ug.^2))/N;