using BoundaryValueDiffEq
using Plots

tstart = 0.0
tstop = 2*pi
tspan = (tstart,tstop)
dt = 0.01
t_lattice = collect(range(tstart,tstop,step = dt))

function simplependulum!(du,u,p,t)
    θ  = u[1]
    dθ = u[2]
    du[1] = dθ
    du[2] = -1*θ
end


function bc!(residual, u, p, t) # u[1] is the beginning of the time span, and u[end] is the ending
    residual[1] = u[1][1]  # the solution at the beginning of the time span should be 0
    residual[2] = u[end][1] # the solution at the end of the time span should be 0
end

bvp = TwoPointBVProblem(simplependulum!, bc!, [0.0,5.0], tspan)
sol = solve(bvp, MIRK4(), dt=dt) # we need to use the MIRK4 solver for TwoPointBVProblem
plot(sol)


function nonConstantInitialGuess(t_lattice)
    θ(t) = sin(t)
    gθ(t) = cos(t)

    vals = Vector{Vector{Float64}}(undef, length(t_lattice))
    for i in 1:length(t_lattice)

        vals[i] = [θ(t_lattice[i]),gθ(t_lattice[i])]
    end

    return vals
end

bvp2 = TwoPointBVProblem(simplependulum!, bc!, nonConstantInitialGuess(t_lattice), tspan)
sol2 = solve(bvp2, MIRK4(), dt=dt) # we need to use the MIRK4 solver for TwoPointBVProblem
plot(sol2)