% One step integration of an ODE by RK2 method, but with the half-step
% slightly advanced (by a factor (1 + epsilon) to stabilize PDEs.
% This is a trick I learned from Dana Longcope.
% CCK 2009-Apr-2
function x1 = rk2plus(ode, x0, t0, h, epsilon)
% f = function handle for ODE,  dx/dt = f(x,t)
% x0 = initial x
% t0 = initial time
% h = timestep (t1 - t0)
% epsilon = timestep bias (<< 1).
% result x1 is x(t1).

xhalf = x0 + (1+epsilon)*(h/2) * ode(x0,t0);
dxdt = ode(xhalf, t0 + (1+epsilon)*h/2);
x1 = x0 + h * dxdt;