NB. Example of simple differential equation.
NB. f=. [: - */  NB. dy/dx = - x*y
NB. Calling params: (x0,y0;x_end;err,htry,hmin)
arg=. ((0 10);1;(1e_8,0.05,1e_8))
kmax=: 100 NB. max number of saved intermediate values
dxsav=: 0.05  NB. minimum save interval
NB. Usage: f sdeint arg -> x_end, y(x_end)

sdeint=: 1 : 0
xx=. x1=. 0{ xy=. xy0=. >0{y
x2=. >1{y
param=. >2{y
eps=. 0{ param
h1=. 1{ param
hmin=. 2{ param
tiny=. 1e_30
maxstp=. 2000
h=. h1 F_sign x2-x1
nok=: nbad=: kount=: 0
xyp=: ''
if. kmax > 0 do. xsav=. xx - 2*dxsav end.
nstp=. 0
while. (nstp=. >:nstp)<: maxstp do. 
 dydx=. u xy
 yscal=. (|}. xy)+(|h*dydx)+tiny
 if. kmax > 0 do.   
  if. (|xx-xsav)>|dxsav do.
   if. kount<(kmax-1) do.
    kount=. >: kount
    xyp=: xyp,xy
    xsav=xx
   end.
  end.
 end.
 if. ((xx+h-x2)*(xx+h-x1))>0 do. h=. x2-xx end.
 out=. u stiff (xy;dydx;h;eps;yscal)
 xx=. 0{ xy=. >0{ out
 hdid=. >1{ out
 hnext=. >2{ out
 if. hdid = h do. nok=: >: nok else. nbad=: >: nbad end.
 if. ((xx-x2)*(x2-x1))>: 0 do.
  if. kmax ~: 0 do. 
   kount=: >: kount
   xyp=: xyp,xy
  end.
  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. Rosenbrock method for stiff equations.    
NB. Numerical Recipes, 2nd ed., p 732.
stiff=: 1 : 0
safety=. 0.9
'grow pgrow shrnk pshrnk errcon'=. 1.5,(%_4),0.5,(%_3),0.1296
maxtry=. 40
xx=. 0{ xysav=. xy=. >0{y
n=. # dysav=. dydx=. >1{y
h=. >2{y
eps=. >3{y
yscal=. >4{y
jout=. jacobn xysav
dfdx=. >0{ jout
dfdy=. >1{ jout
jtry=. 0
while. (jtry=. >: jtry) <: maxtry do.
 ai=. %. ((2%h)* = i. n) - dfdy  NB. inverse of matrix
 g1=. ai +/ .* dysav + h*dfdx%2
 xy=. xysav + (h,(2*g1))
 dydx=. u xy
 g2=. ai +/ .* dydx + (h*_1.5*dfdx) + _8*g1%h
 xy=. xysav + ((0.6*h),((48*g1%25)+(6*g2%25)))
 dydx=. u xy
 g3=. dydx + (h*121*dfdx%50) + ((372*g1%25)+(12*g2%5))%h 
 g3=. ai +/ .* g3 
 g4=. dydx+(h*29*dfdx%250)+((_112*g1%125)+(_54*g2%125)+(_2*g3%5))%h
 g4=. ai +/ .* g4
 xy=. xysav +(h,((19*g1%9)+(g2%2)+(25*g3%108)+(125*g4%108)))
 err=. (17*g1%54)+(7*g2%36)+(125*g4%108)
 if. (0{xy) = (0{xysav) do. 
   'Step too small!' 
   return. 
 end.
 errmax=.(%eps)* >./ | err%yscal
 if. errmax <: 1 do. 
  if. errmax > errcon do.
   hnext=. safety*h*(errmax^pgrow)
  else.
   hnext=. grow*h
  end.
  z=. (xy;h;hnext)
  return.
 else.
  hnext=. safety*h*(errmax^pshrnk)
  h=.((|hnext)>. (shrnk*|h)) F_sign h
 end.
end.
'Exceeded MAXTRY in stiff.'
)

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

NB.  A simple stiff system:

NB. A=: 2 2 $ 998 1998 _999 _1999
NB. L=: [ +/ .* }.@]
NB. f =: A&L

NB. jacobn=: 3 : 0
NB. xy=. y          NB. not used - const. coeff.
NB. dfdx=. 2 # 0
NB. dfdy=. 2 2 $ 998 1998 _999 _1999
NB. dfdx;dfdy
NB. )

NB. Example of stiff set of 3 equations:
G=: 3 : 0
'x y1 y2 y3'=. y
dy1=. (_0.013*y1)+(_1000*y1*y3)
dy2=. _2500*y2*y3
dy3=. (_0.013*y1)+(_1000*y1*y3)+(_2500*y2*y3)
dy1,dy2,dy3
)

NB. The Jacobian for the above system.
jacobn=: 3 : 0
dfdx=. 3 # 0
'x y1 y2 y3'=. y
dfdy=. 3 3 $ 0
dfdy=. ((_0.013-1000*y3),0,(_1000*y1)) 0}dfdy
dfdy=. (0,(_2500*y3),(_2500*y2)) 1}dfdy
a1=. _0.013-1000*y3
a2=. _2500*y3
a3=. (_1000*y1)+(_2500*y2)
dfdy=. (a1,a2,a3) 2}dfdy
dfdx;dfdy
)

NB. Integrate equations from 0 to 5:
   G sdeint ((0 1 1 0);5;(1e_4 2.9e_4 1e_10))
5 0.954056 1.04594 _3.47517e_6
   nok
38
   nbad
2

NB. Try same with Runge-Kutta (runge_kutta.ijs):
   G odeint ((0 1 1 0);5;(1e_4 2.9e_4 1e_10))
5 0.954056 1.04594 _3.47519e_6
   nok
3717
   nbad
1010