# Solving Boundary Value Problems with Shooting() and DifferentialEquations.jl

Hey, Im trying to solve a Boundary Value Problems with Shooting() and DifferentialEquations.jl.

My problem is formulated as this:

This problem is an IVP with known initial values and has also a known analytical solution.
The problem is solved as an ODE:

``````using DifferentialEquations

const μ = 0.04                                                   # equity premium
const σ = 0.13                                                   # volatility of equity return (std)
const γ = 5.0                                                    # risk aversion parameter
const ρ = 0.03                                                   # discount rate
const r = 0.02                                                   # risk-free rate

t_min = 1.0                                                      # start
t_max = 50.0                                                     # end
tspan = (t_min,t_max)                                            # fund value interval
step_size = 1                                                    # saved solutin every step of size

η = (1/γ)ρ+(1-1/γ)*(r+(1/2)*(μ^2/(γ*σ^2)))                       # calculate eta

# Defining the IVP and solving it.
x₀ = η^(-γ)*(t_min^(1-γ)/(1-γ))                                                                   # initial value for the value function
y₀ = η^(-γ)*t_min^(-γ)                                                                            # initial value for the first derivative of the value function
u₀ = [x₀, y₀]                                                                                     # initial state vector

function nohabits!(du,u,p,t)
V = u[1]
V_a = u[2]
du[1] = V_a                                                                                  # the first derivative of the value (dx) function equals y
du[2] = (0.5*V_a^2*(μ/σ)^2)/((γ/(1-γ))*V_a^(1-1/γ)-ρ*V+r*t*V_a)                              # equation 12a from the notes
end

prob = ODEProblem(nohabits!,u₀,tspan,maxiters=10e5)                                              # define the problem
@time sol = solve(prob,reltol=1e-12,abstol=1e-12, saveat=step_size)                              # solve the problem
``````

The solution time is: 35.300468 seconds (86.01 M allocations: 4.058 GiB, 6.43% gc time) for the first run. And for the 2nd run: 0.000710 seconds (881 allocations: 75.578 KiB).

If the tolerance level is increased to reltol=1e-8, abstol=1e-8, the error becomes present.
And at a tolerance level of reltol=1e-9,abstol=1e-9 it becomes insignificant.

My true problem has the same form as the previous, but the initial conditions are not known.
Hence, I want to solve the previous problem as a BVP, using Shooting() to find the initial conditions before attempting to solve my true problem.

This problem is formulated as this:

``````# Building  and solve the BVP
function nohabits2!(du,u,p,t)
du[1] = u[2]
du[2] = (0.5*u[2]^2*μ^2/σ^2)/((γ/(1-γ))*max(u[2],0)^(1-1/γ)-ρ*u[1]+r*t*u[2])
end

off_set = 0
x_guess = x₀ + off_set                      # My guess of V_init
V_a_end = sol.u[end][2]                     # Last value of V_a in the known solution
u₀_2 = [x_guess, y₀]                        # The initial guess

function bc!(residual,u,p,t)
residual[1] = u[end][2] - V_a_end
end

bvp3 = BVProblem(nohabits2!, bc!, u₀_2, tspan, maxiters=1e5)
@time sol3 = solve(bvp3, Shooting(Tsit5()), isoutofdomain=(u,p,t) -> u[2]<0, saveat=step_size)
``````

And solved by guessing the correct initial conditions.
This problems solves in: 10.299134 seconds (8.04 M allocations: 378.491 MiB, 1.45% gc time)
And, retcode: MaxIters.

By increasing maxiters=1e7
This challenge diapers, and the problem solves quickly but with an high error:
10.142641 seconds (3.72 M allocations: 178.774 MiB, 0.46% gc time)
retcode: Success

When trying to reduce the error by reducing the tolerance level and max iterations: maxiters=1e8 and reltol=1e-5,abstol=1e-5. The solution times become: 199.300782 seconds (3.72 M allocations: 178.776 MiB, 0.02% gc time)

A tolerance level at reltol=1e-5,abstol=1e-5 yields a high error. And from the first problem we know that the tolerance level should be reltol=1e-9,abstol=1e-9.

In the search for a proper solution I have tried several other algorithms, such as: `AutoTsit5(Rosenbrock23()), AutoVern7(Rodas5()), Vern7(), Vern9() ` etc.
with no success.

So the question is:
Why is the BVP so hard to solve when the initial guess is correct?