RungeKutta4.hs

-- https://www.wolframalpha.com/input/?i=y%27%3D-2*x%2B3*y-5

(t0,y0) = (0,1)
h = 0.01

type Function = Double -> Double

yt :: Function
yt t = (y0 - 17/9) * (exp $ 3*t) + 2*t/3 + 17/9

type Function2 = Double -> Double -> Double

yt1 :: Function2
yt1 t y = -2*t + 3*y - 5

rk4 :: Function2 -> Double -> (Double,Double) -> (Double,Double)
rk4 f h (tn,yn) =
 (tn1,yn1)
 where
  k1 = f tn yn
  k2 = f (tn + h/2) (yn + h*k1/2)
  k3 = f (tn + h/2) (yn + h*k2/2)
  k4 = f (tn + h) (yn + h*k3)
  tn1 = tn + h
  yn1 = yn + (k1 + 2*k2 + 2*k3 + k4) * h/6

printIteration (tn,yn) =
 putStrLn string >>
 return (tn1,yn1)
 where
  string =
   "t = " ++ show tn ++ ", RK4 y = " ++ show yn ++ ", y = " ++ show y ++
   ", error = " ++ show error
  (tn1,yn1) = rk4 yt1 h (tn,yn)
  y = yt tn
  error = abs $ y - yn

printIterations (tn,yn) finalT =
 printIteration (tn,yn) >>= \(tn1,yn1) ->
 if tn1 > finalT then
  return ()
 else
  printIterations (tn1,yn1) finalT

main =
 printIterations (t0,y0) 5