struct Problem
    # u'' = -g*u + s
    g
    s
    tstops
end

mutable struct Cache
    u
    uprev
    t
    tprev
    tnext
    h # stepwidth
    step # current step
    sol # array of solutions
end


function solve(prob,u0,du0)

    cache = Cache(undef,undef,undef,undef,undef,undef,undef,undef)
    initialize(prob,cache,u0,du0)

    for t in 3:length(prob.tstops)
        perform_step(prob,cache)
    end

    return cache

end

function initialize(prob,cache,u0,du0)
    cache.sol = Array{Float64}(undef,length(prob.tstops))
    cache.tprev = prob.tstops[1]
    cache.t = prob.tstops[2]
    cache.tnext = prob.tstops[3]
    cache.h = cache.t-cache.tprev
    cache.u = u0 +
              du0*cache.h +
              cache.h^2 / 2 * (-prob.g(cache.tprev)*u0 + prob.s(cache.tprev)) +
              cache.h^3 /6 * (-prob.g(cache.tprev) )

    cache.uprev = u0
    cache.sol[1] = cache.uprev
    cache.sol[2] = cache.u
    cache.step = 3
end

function perform_step(prob,cache)
    u = cache.u
    uprev = cache.uprev

    t = cache.t
    tprev = cache.tprev
    tnext = cache.tnext

    g = prob.g(t)
    gnext = prob.g(tnext)
    gprev = prob.g(tprev)

    s = prob.s(t)
    snext = prob.s(tnext)
    sprev = prob.s(tprev)

    h = cache.h

    unext = (2*u*(1-(5*h^2*g)/12) -
    uprev*(1+(h^2*gprev)/12) +
    h^2/12*(snext+10*s+sprev))/
    (1+h^2/12*gnext)

    cache.uprev = cache.u
    cache.u = unext
    cache.sol[cache.step]=unext


    cache.step = cache.step+1
    step = cache.step

    if step >= length(prob.tstops)
        return
    end
    cache.tprev = t
    cache.t = tnext
    cache.tnext = prob.tstops[step]
end


function main(g,s)
	#u'' = -gu + s
	#
    #g = 1
    #s = 8
    tstops = 0:0.0001:10
    prob = Problem(x->g,x->s,tstops)
    @time sol = solve(prob,0,0)
    unicodeplots()
    y = (broadcast(cos,prob.tstops.*sqrt(g)).-1) .*(-s/g)
    plot(prob.tstops,sol.sol)
    plot!(prob.tstops,y)
    xlims!(prob.tstops[1],prob.tstops[end])
    gui()
end

main(rand(2)...)