Hello,

I am trying to replicate the solution of an optimal control problem (Redirecting).

The canonical system is given by

where x is the state variable, \lambda the costate, and u the control.

The maximizer is given by u=\frac{x}{2}(p-\lambda).

The boundary conditions are the initial condition x(0) =x_0 and the transversality condition \lim_{t\rightarrow \infty} e^{-rt}\lambda(t)=0.

This system has a hyperbolic saddle point \hat y and I am trying to solve for the stable path from the initial state to the steady state.

Let

and denote by J(\hat y) the Jacobian matrix at the the steady state.

Let E^s(\hat y) be the eigenspace associated with the negative eigenvalues, which is tangent to the stable manifold at \hat y.

Then, for T large enough, the transversality condition can be approximated by F^\prime(y(T)-\hat y) =0, where F^\prime is a matrix spanning the orthogonal complement to the stable eigenspace.

My try of implementing this so far:

```
using DifferentialEquations
using Optim
using ForwardDiff
using LinearAlgebra
using Plots
r = 0.1
τ = 0.25
c = 0.375
a = 7.5
pr = 0.1
par = [r, τ, c, a, pr]
T = 500
G(X) = X * (1 - X / a)
P(X) = (1 / (c + τ)) * (X^2 / (1 + X^2))
```

Getting the Jacobian matrix and F

```
function f(u)
fish = u[1] / 2 * (pr - u[2])
du1 = G(u[1]) - P(u[1]) - u[1] * fish
du2 = u[2] * (r - 1 + 2 * u[1] / a + 2 * u[1] / ((c + τ) * (1 + u[1]^2)^2)) - fish * (pr - u[2])
return [du1, du2]
end
Df(u) = ForwardDiff.jacobian(f, u)
function get_steady_state(par)
r, τ, c, a, pr = par
fish(x, barλ) = x * (pr - barλ) / 2
eqx(x, λ) = G(x) - P(x) - x * fish(x, λ)
eqλ(x, λ) = λ * (r - 1 + 2 * x / a + 2 * x / ((c + τ) * (1 + x^2)^2)) - fish(x, λ) * (pr - λ)
eqs(x, λ) = [eqx(x, λ), eqλ(x, λ)]
res = Optim.optimize(x -> sum(abs.(eqs(x[1], x[2]))), [3, 0.5], Optim.Options(g_tol=1e-15))
x = res.minimizer[1]
λ = res.minimizer[2]
return x, λ
end
barx, barλ = get_steady_state(par)
J = Df([barx, barλ])
ev = eigvals(J)
stable_evecs = eigvecs(J)[:, ev.<0]
F = nullspace(stable_evecs')
```

Boundary conditions

```
x0 = barx + 0.1
function bc1!(residuala, sola, p)
par, x0, barx, barλ = p
residuala[1] = sola[1][1] - x0
end
function bc2!(residualb, solb, p)
par, x0, barx, barλ = p
residualb[1] = (F'*(solb[end].-[barx, barλ]))[1]
end
```

Here I also ran into a problem when trying to call `sola(0.0)[1]`

as in the tutorial (Boundary Value Problems · DifferentialEquations.jl) it returns an error

```
ERROR: MethodError: objects of type Vector{Float64} are not callable
Use square brackets [] for indexing an Array.
Stacktrace:
[1] bc1!(residuala::SubArray{…}, sola::Vector{…}, p::Tuple{…})
```

Solving

```
strt = [barx, barλ]
prob = TwoPointBVProblem(f!, (bc1!, bc2!), strt, (0.0, 500.0), (par, x0, barx, barλ); bcresid_prototype=(zeros(1), zeros(1)))
@time sol = solve(prob, MIRK4(), dt=0.01)
```

Which yields

where for some reason the solution diverges at the end, violating the boundary condition.

Help would be appreciated!