Just a quick update: ApproxFun works really well for this problem. Here’s the code I used to solve it. It’s very similar to @ptoche’s response using DifferentialEquations.
using ApproxFun
using Plots
maxt = 100 # Time periods
x=Fun(identity, 0..maxt) # Let x() be an empty Chebychev fit on the range
# [0, 100]
#Parameters
ρ = 0.10
δ = 0.10
α = 0.5
A = 0.4
σ = 0.5
p = [ρ, δ, α, A, σ]
Y(K) = A * K ^ α - δ * K # output
#Steady-States
K_ss = ((δ + ρ) / (A * α)) ^ (1 / (α - 1))
C_ss = (A * K_ss ^ α - δ * K_ss)
#Starting from above steady-state
ramseyhigh = (K,C) -> [K' - (A * K^α - δ * K - C );
C' + C * σ * (ρ + δ - α * A * K ^ (α - 1));
K(0) - K_ss * 10.0;
K(maxt) - K_ss * 1.00001]
#Starting from below steady-state
ramseylow = (K,C) -> [K' - (A * K^α - δ * K - C );
C' + C * σ * (ρ + δ - α * A * K ^ (α - 1));
K(0) - K_ss * 0.01;
K(maxt) - K_ss * 0.99999]
# Choosing a starting function for capital
# and consumption.
# The approximation would probably be
# faster If I had a better starting function.
K10 = one(x) * K_ss
C10 = one(x) * C_ss
# Run the function approximation
Kh,Ch = newton(ramseyhigh, [K10,C10], maxiterations=50, tolerance = tolerance=1E-10)
Kl, Cl = newton(ramseylow, [K10,C10], maxiterations=50, tolerance = tolerance=1E-10)
#Recreating the plot from ptoche's post:
plt = plot(title="Policy Function in Ramsey-Cass-Koopmans Model", xlims=(0,2), ylims=(0,0.6), xlabel="K(t)", ylabel="C(t)", legend=false)
plot!(plt, x -> Y(x), 0, 2, color="black", linewidth=3, legend=false)
vline!(plt, [K_ss], color="black", linewidth=3, legend=false)
plot!(plt, Kl, Cl, xlims=(0,2), ylims=(0,0.6), color=:blue, linewidth=2)
scatter!((K_ss,C_ss), label=["Stationary State" "."], color=:red)
plot!(plt, Kh, Ch, xlims=(0,2), ylims=(0,0.6), color=:blue, linewidth=2)
savefig(plt, "ramsey-approx_fun.png")
