# L-31 MCS 471 Wed 2 Nov 2022 : rk45stiff.jl

# Using RK45 of the scipy.integrate module on a stiff problem.

using SciPy
using Printf

println("Solving a stiff problem with SciPy.integrate.RK45")

"""
    rhs(t, y) 

Defines the right hand side of the ODE.
"""
function rhs(t, y)
    result = zeros(length(y))
    result[1] = 10*(1.0 - y[1])
    return result
end

sol = SciPy.integrate.RK45(rhs, 0.0, (0.5, ), 100.0, rtol=1.0e-8)

println("step :     h    :    t     :           y            :  error")

stepcnt = 0
while sol.status != "finished"
    sol.step()
    global stepcnt = stepcnt + 1
    scnt = @sprintf("%4d", stepcnt)
    step = @sprintf("%.2e", sol.step_size)
    tval = @sprintf("%.2e", sol.t)
    soly = @sprintf("%.16e", sol.y[1])
    exact = 1 - exp(-10*sol.t)/2
    err = abs(sol.y[1] - exact)
    serr = @sprintf("%.2e", err)
    println("$scnt : $tval : $step : $soly : $serr")
end