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 |