Hi all,

I am interested in using the `BoundaryValueDiffEq`

part of `DifferentialEquations`

to solve economic growth problems in continuous time. The problem is, I can’t get a relatively simple example to work!

The example I am trying to solve is the Cass-Koopmans-Ramsey model from Aghion and Howitt’s *Endogenous Growth Theory*. It is a nonlinear system of differential equations

With an initial capital stock K(0) = K_0 and an end-of-time condition

\alpha, \delta, and \rho are parameters in the model and represent the productivity of capital, the depreciation rate of capital, and the intertemporal preference of a utility maximizer.

The steady-state for this problem is easy to solve for analytically, but I’m interested in what happens along the transitional path, hence, the need for numerical analysis.

Recognizing that `BoundaryValueDiffEq`

will not appreciate a boundary condition at t = \infty, I made the following transformation:

which transforms my system to

and my terminal condition to

The terminal condition is still somewhat convoluted, so I can exploit what I know about the steady-state to replace it. I know that the system should reach the steady-state at \tau = 1, so, after solving for the steady-state, I replace my terminal condition with

where

Now that my system is defined to my satisfaction, I turn to the code:

```
using DifferentialEquations
using Plots
#Parameters
ρ = 0.04
δ = 0.10
α = 0.66
p = [ρ, δ, α]
#Steady-States
K_ss = ((δ + ρ) / α) ^ (1 / (α - 1))
C_ss = (K_ss ^ α - δ * K_ss)
#Initial Capital
K0 = 90.0
#Defining the system
function end_sav!(du, u, p, t)
ρ = p[1]
δ = p[2]
α = p[3]
K = u[1]
C = u[2]
du[1] = (K^α - δ * K - C ) / (1 - t) ^ 2
du[2] = - C * (ρ + δ - α * K ^ (α - 1)) / (1 - t) ^ 2
end
# Define the boundary conditions
function bc!(residual, u, t)
residual[1] = u[1][1] - K0
residual[2] = u[2][end] - C_ss
end
# Prevent the solver from trying
# negative values of capital and consumption
function g(u, resid)
resid[1] = u[1] - 0.001
resid[2] = u[2] - 0.001
end
end_sav_prob = BVProblem(end_sav!, bc!, [0.90, C_ss], (0.0, 1.), p; callback = GeneralDomain(g))
sol = solve(end_sav_prob)
```

I’ve tried different algorithms, different starting values, different time scales, varying the parameters, etc. Nothing seems to work. The most common issue is instability. Any suggestions on what’s going on? Additionally, I can’t be the first person to try to use Julia to solve these kinds of problems. Does anyone know of any resources/packages for working with continuous-time models in Julia? QuantEcon has some really good stuff, but as far as I can tell they stick to discrete-time methods.