module ODE
let isBetween (a1, a2) x =
if a1 < a2 then a1 <= x && x<=a2
else a2 <= x && x<=a1
let SolveEuler (f:float*float->float) (h:float) (init:float*float) (T:float) =
let t0, _ = init
Seq.unfold (fun ty ->
let t_, y_ = fst ty, snd ty
if not (isBetween(t0, T) t_) then None
else
Some (ty, (t_ + h, y_ + h*(f ty)))
) init
let SolveRungeKutta (f:float*float->float) (h:float) (init:float*float) (T:float) =
let t0, _ = init
Seq.unfold (fun ty ->
let t_, y_ = fst ty, snd ty
if not (isBetween(t0, T) t_) then None
else
let k1 = f ty
let k2 = f (0.5*h+t_, y_+0.5*h*k1)
let k3 = f (0.5*h+t_, y_+0.5*h*k2)
let k4 = f ( h+t_, y_+h*k3)
let t,y = t_+h, y_ + h*(k1+2.0*k2+2.0*k3+k4)/6.0
Some (ty, (t, y))
) init |