NB. Runge Kutta differential equation integrator.
NB. Example of simple differential equation:
NB. dy/dx = - x*y  we write as ff=: 3 : '-(0{y)*(1{y)'
NB. Calling params: x0,y0;x_end;err,htry,hmin
NB. e.g.  arg=. (0 10);1;(1e_8,0.05,1e_8)
NB. Usage: ff odeint arg --> x_end, y(x_end)

odeint=: 1 : 0
'xy0 x2 param'=. y
'eps h1 hmin'=. param
x=. x1=. 0{ xy=. xy0
tiny=. 1e_30
maxstp=. 10000
h=. h1 F_sign x2-x1
nstp=. nok=: nbad=: kount=: 0
while. (nstp=. >:nstp)<: maxstp do. 
 dydx=. u xy
 yscal=. (|}. xy)+(|h*dydx)+tiny
 if. ((x+h-x2)*(x+h-x1))>0 do. h=. x2-x end.
 'xy hdid hnext'=. u rkqs xy;dydx;h;eps;yscal
 x=. 0{ xy
 if. hdid = h do. nok=: >: nok else. nbad=: >: nbad end.
 if. ((x-x2)*(x2-x1))>: 0 do.
  xy0=. xy
  return.
 end.
 if. (|hnext)<hmin do.
  'Stepsize smaller than minimum in odeint.'
  return.
 end.
 h=. hnext
end.
'Too many steps in odeint.'
)

NB. Cash-Karp Runge-Kutta quality control step
NB. Numerical Recipes, 2nd ed., p 713.
rkqs=: 1 : 0
'safety pgrow pshrink errcon'=. 0.9 _0.2 _0.25 1.89e_4
'xy dydx h eps yscal'=. y
whilst. 1 do.
 'xynew yerr'=. u rkck xy;dydx;h
 errmax=: >./ | yerr%yscal*eps
  if. errmax > 1 do. 
  htemp=. safety*h*errmax^pshrink
  h=.((|htemp)>. (0.1*|h)) F_sign h
  xnew=. (0{xy)+h
   if. xnew = (0{xy) do. 
    'Step too small!' 
    return. 
   end.
  continue.
  else. 
   if. errmax > errcon do.
    hnext=. safety*h*errmax^pgrow
   else.
    hnext=. 5*h
   end.
   xynew;h;hnext
  end.
  return.
end.
)

NB. Cash-Karp Runge-Kutta (5th order) adverb. Usage:        
NB. f rkck xy0;dydx0;h --> xy1=[x+h,y(x+h)];error
rkck=: 1 : 0
'xy0 dydx0 h'=. y
k1=. h * dydx0
k2=. h * u xy0+(h,k1)%5
k3=. h * u xy0+((0.3*h),((3*k1%40)+(9*k2%40)))
dy=. ((3*k1)+(_9*k2)+(12*k3))%10
k4=. h * u xy0+((0.6*h),dy)
dy=. (_11*k1%54)+(2.5*k2)+(_70*k3%27)+(35*k4%27)
k5=. h * u xy0+(h,dy)
dy=. (1631*k1%55296)+(175*k2%512)+(575*k3%13824)+(44275*k4%110592)+(253*k5%4096)
k6=. h * u xy0+((0.875*h),dy)
dy=. (37*k1%378)+(250*k3%621)+(125*k4%594)+(512*k6%1771) NB. 5th order step
xy=. xy0+(h,dy)
dy4=. (2825*k1%27648)+(18575*k3%48384)+(13525*k4%55296)+(277*k5%14336)+(k6%4)
xy;(dy-dy4)  NB. dy4 is embedded 4th order step; (dy-dy4) is the error estimate 
)

NB. Model of FORTRAN "SIGN(A,B)" function:
F_sign=: ([: | [)* 1: -[: +: ]< 0: