The infinite-horizon continuous-time problems with a saddle-point structure could be solved using the time reversal technique outlined here. At some point some more general technique should be put together though. 
Edit Adding some code for the simple deterministic saddle-path. As it turns out, the optimal trajectory features a constant level of investment (if I copied the equations above correctly); for this reason I have drawn the optimal trajectory with dashes, so we can see the black nullcline below it. Solving this problem is just like solving the Ramsey-Cass-Koopmans model (see the link above). Both models should be part of a general setup, as discussed below by Chris and Jesse.
using DifferentialEquations
using Plots
# Parameters
p = Dict(:δ => 0.1, :ρ => 0.05, :z => 0.3)
"""firm investment problem
K: state variable, with given initial value K0
I: control variable, free to jump"""
function inv!(du,u,p::Dict,t)
δ, ρ, z = p[:δ], p[:ρ], p[:z]
K, I = u
du[1] = I - δ * K
du[2] = (ρ + δ - z) + (ρ + δ)*I
nothing
end
"""compute the stationary state"""
function ss(p::Dict)
δ, ρ, z = p[:δ], p[:ρ], p[:z]
Iss = z / (ρ+δ) - 1
Kss = Iss/δ
return Dict(:Kss => Kss, :Iss => Iss)
end
# The stationary state
SS = ss(p) # returns 10.0 and 1.0 for reference parameter values
Kss, Iss = ss(p)[:Kss], ss(p)[:Iss]
# event handler to prevent variables going out of bounds
function set_callback()
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*Iss - 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
cb = set_callback()
# Run time "backwards" | start below the stationary state
tspan = (0.0, -100.0)
tweak = 1.01 # I0 starts above Iss
u0 = [0.999999 * Kss; 0.999999 * tweak * Iss]
ode = ODEProblem(inv!, u0, tspan, p)
sol1 = solve(ode, abstol=1e-8, reltol=1e-8, callback=cb)
# Run time "backwards" | start above the stationary state
tspan = (0.0, -100.0)
tweak = 0.99 # I0 starts below Iss
u0 = [1.00001 * Kss; 1.00001 * tweak * Iss]
ode = ODEProblem(inv!, u0, tspan, p)
sol2 = solve(ode, abstol=1e-8, reltol=1e-8, callback=cb)
# initialize plot
plt = plot(title="Policy Function in the Neoclassical Investment Model", xlims=(0,2*Kss), ylims=(0,2*Iss), xlabel="K(t)", ylabel="I(t)", legend=false)
# plot nullclines
plot!(plt, x-> p[:δ] * x, 0, 2*Kss, color=:black, linewidth=3, legend=false)
plot!(plt, x-> Iss, 0, 2*Kss, color=:black, linewidth=3, legend=false)
plot!(plt, sol1, xflip=false, vars=(1,2), xlims=(0,2*Kss), ylims=(0,2*Iss), color=:blue, linewidth=2, linestyle=:dash)
plot!(plt, sol2, xflip=false, vars=(1,2), xlims=(0,2*Kss), ylims=(0,2*Iss), color=:blue, linewidth=2, linestyle=:dash)
# add decoration
scatter!((Kss,Iss), label=["Stationary State" "."], color=:red)
savefig(plt, "investment-saddle.png")
# make the streamplot
using CairoMakie
"""firm investment model equations as points in (K,I) space"""
function invEq(K, I, p)
δ, ρ, z = p[:δ], p[:ρ], p[:z]
dK = I - δ * K
dI = (ρ + δ - z) + (ρ + δ)*I
return Point(dK, dI)
end
let
fig = Figure(fontsize = 20)
ax = Axis(fig, xlabel="K", ylabel="I", title="Streamplot of Neoclassical Investment Model", backgroundcolor=:black)
streamplot!(ax, (x,y)->invEq(x,y,p), 0..2*Kss, 0..2*Iss, colormap = Reverse(:plasma),
gridsize= (32,32), arrow_size = 10)
fig[1, 1] = ax
display(fig)
save("investment-streamplot-black.png", fig)
end
