The following two-point boundary value problem is difficult to solve because the stationary state has a saddle structure, with one positive and one negative eigenvalue: only one path leads to the stationary state — all other paths are unstable. Can the `BVProblem`

options be tweaked to achieve convergence?

```
using DifferentialEquations
using Plots
"""Parameters (including the stationary state values)"""
p = Dict(:σ => 0.5, :K0 => 0.2, :C0 => 0.1, :Kss => 1, :Css => 0.3)
"""The Production function"""
F(x) = 0.4 * x^(0.5) - 0.1 * x
"""The Ramsey-Cass-Koopmans Dynamic system"""
function sys!(du,u,p::Dict,t)
σ = p[:σ]
K, C = u
du[1] = 0.4 * K^(0.5) - 0.1 * K - C
du[2] = (0.5 * 0.4 * K^(-0.5) - 0.1 - 0.1) * σ * C
end
"""boundary condition
The state variable K has given initial value K0.
The control variable C has known stationary-state value Css."""
function bc!(res,sol,p::Dict,t)
K0, Css = p[:K0], p[:Css]
res[1] = sol[1][1] - K0
res[2] = sol[end][1] - Css
end
# event handler to prevent variables going out of bounds
function set_callback(p::Dict)
Kss, Css = p[:Kss], p[:Css]
c1(u, t, integrator) = u[1] # stay in positive range
c2(u, t, integrator) = u[2] # stay in positive range
c3(u, t, integrator) = 2*Kss - u[1] # don't stray too far from stationary state
c4(u, t, integrator) = 2*Css - u[2] # don't stray too far from stationary state
function terminate_affect!(integrator)
terminate!(integrator)
end
# set up affect upon callback
cb1 = ContinuousCallback(c1, terminate_affect!)
cb2 = ContinuousCallback(c2, terminate_affect!)
cb3 = ContinuousCallback(c3, terminate_affect!)
cb4 = ContinuousCallback(c4, terminate_affect!)
cb = CallbackSet(cb1, cb2, cb3, cb4)
return cb
end
# SOLVE AS TWO-POINT BVP
const u0 = [p[:K0]; p[:C0]]
const tspan=(0, 5)
bvp = BVProblem(sys!, bc!, u0, tspan, p)
cb = set_callback(p)
# attempt to solve bvp with shooting method
sol = solve(bvp, callback=cb, Shooting(Vern7()))
# PHASE PLOT
plt = plot(xlims=(0,2*p[:Kss]), ylims=(0,2*p[:Css]), xlabel="K(t)", ylabel="C(t)",
extra_kwargs=Dict(:subplot=>Dict(:legend_hfactor=>1.5)))
# plot the nullclines
plot!(plt, F, 0, 2*p[:Kss], color=:black, linewidth=2, label="")
vline!(plt, [p[:Kss]], color=:black, linewidth=2, label="")
# mark the position of the stationary state
scatter!(plt, (p[:Kss],p[:Css]), label=["Stationary State" "."], color=:red)
# mark the position of the initial guess u0
scatter!(plt, (u0...,), label=["Initial Guess" "."], color=:green)
# add trajectory to plot
plot!(plt, sol, vars=(1,2), xlims=(0,2*p[:Kss]), ylims=(0,2*p[:Css]), color=:blue, linewidth=3, label="converging path", thickness_scaling=1)
display(plt)
```

This is a follow-up on solving-boundary-value-differential-equation-problems-for-economics. I decided to start a new topic since the problem was solved there via other methods.

Can the `BVP`

solver find the stable path? This is what the path that converges to the saddle looks like.

And this is what my attempts to use `BVProblem`

yielded: