# Plotting Phase Diagram

I have the following problem that I'm sure Mathematica can handle, but it's not working for me!

In the following code, I'm trying to replicate the Ramsey Model Phase Diagram.

In fact, I want now to be able to dram the stable arm of the system. My attempted is using NDSolve, since I cannot do better. The idea is given a $z_0$ in the regions where exist a stable arm (the ones where the blue line cuts). I tried the framework of BV problems and add up some 1% range; I'm rather sure Mathematica can do better than this.

Could anybody suggest me a way to solve it, because NDSolve refuses to work with this coding, and I really don't know how to do it.

α = 0.3;
ρ = 0.02;
θ = 2;
η = 0.5;
δ = 0.05;
A = 1.0;

T = 600;
err = 0.01;
{zM, cM} = {250000, 1000};

Dz = A α^α η^η (α+η)^-(α+η) z^(α+η)-δ z- c;
Dc = θ^-1 (A α^α η^η (α+η)^(1-(α+η))  z^(α+η-1) - ρ-δ);
{zT, cT} = NSolve[{Dc == 0, Dz == 0}, {z, c}][[1, All, 2]]

Plot1 = ContourPlot[{Dz==0,Dc==0},
{z, 0, zM}, {c, 0, cM},
PlotTheme -> "Scientific",
PlotLegends -> {"Subscript[Overscript[z, .], t]=0", "Subscript[Overscript[c, .], t]=0"},
Frame -> False, Axes -> True, AxesLabel -> {"z", "c"},
PlotRange -> {{0, zM}, {0, cM}}, ImageSize -> Large
];

z0 = 10000;

sol=NDSolve[
{z'[t] == A α^α η^η (α+η)^-(α+η) z[t]^(α+η)- δ z[t]- c[t],
c'[t] == c[t]/θ (A α^α η^η (α+η)^(1-(α+η)) z[t]^(α+η-1) - ρ-δ),
Abs[z[T]-zT]<err*zM,
Abs[c[T]-cT]<err*cM
},
{z, c}, {t, 0, T},
Method -> {"Shooting", "StartingInitialConditions" -> {z[0] == z0}}
];

Plot2 = ParametricPlot[Evaluate[{z[t], c[t]} /.sol], {t,0,T},
PlotRange -> {{0, zM}, {0, cM}}, ImageSize -> Large
];

Show[{Plot1, Plot2}]

• This is pretty much a code dump. Could you at least remove the useless formatting in your plots, reduce the code down to a minimal example, indicate what happens when you try to run NDSolve, whether you get an error etc. Why should we do all that for you? – MarcoB May 10 '16 at 18:09
• What do you mean by "stable arm of the system"? – Chris K May 10 '16 at 22:55
• You can ameliorate my solution by adding to ContourPlot PlotPoints -> 50. – cyrille.piatecki May 11 '16 at 7:23
• When I export in eps that graphic seems perfect. The problem is only with png – cyrille.piatecki May 11 '16 at 7:24

I am badly placed to critic you question since I use to be as criptic as you. First only an economist can understand your objective since this representation is not so standard as you think. I have tried a day long to change the parameters but it is nearly impossible to change the position of c=0. Using the free add-on of Gianluca Gorni here is the best I have been able to do --- I have changed the font to times since I guess you are typesetting with Tex.

Here is the code.

Needs["CurvesGraphics6"];
ClearAll[α, ρ, θ, η, δ, A, z, c, t];
α = 0.3;
ρ = 0.001;
θ = 2.;
η = 0.09;
δ = 0.05;
A = 1;
Φ =
A α^α η^η (α + η)^-(α + \
η);
Ω =
A α^α η^η (α + η)^(
1 - (α + η));
{zM, cM} = {250., 1.}
F[z_, c_] := Φ z^(α + η) - δ z - c
G[z_, c_] :=
1/θ (Ω z^((α + η) -
1)) - (ρ + δ)
NSolve[Φ z^(α + η) - δ z == 0, z]
NSolve[1/θ (Ω z^((α + η) -
1)) - (ρ + δ) == 0, z]
sol1 = NDSolve[( {
{z'[t] == Φ z[t]^(α + η) - δ z[
t] - c[t]},
{c'[t] ==
1/θ (Ω z[t]^((α + η) -
1)) - (ρ + δ)},
{z[0] == 5},
{c[0] == 1.5}
} ), {z, c}, {t, 0, 800}, PlotSolution -> True,
ColorFunction -> "DarkRainbow", PlotPoints -> 50];
sol2 = NDSolve[( {
{z'[t] == Φ z[t]^(α + η) - δ z[
t] - c[t]},
{c'[t] ==
1/θ (Ω z[t]^((α + η) -
1)) - (ρ + δ)},
{z[0] == 15.6342},
{c[0] == 2}
} ), {z, c}, {t, 0, 100}, PlotSolution -> True,
ColorFunction -> "DarkRainbow"];
sol3 = NDSolve[( {
{z'[t] == Φ z[t]^(α + η) - δ z[
t] - c[t]},
{c'[t] ==
1/θ (Ω z[t]^((α + η) -
1)) - (ρ + δ)},
{z[0] == 20},
{c[0] == 1.2}
} ), {z, c}, {t, 0, 30}, PlotSolution -> True,
ColorFunction -> "DarkRainbow"];
sol4 = NDSolve[( {
{z'[t] == Φ z[t]^(α + η) - δ z[
t] - c[t]},
{c'[t] ==
1/θ (Ω z[t]^((α + η) -
1)) - (ρ + δ)},
{z[0] == 0.429975},
{c[0] == .2}
} ), {z, c}, {t, 0, 80}, PlotSolution -> True,
ColorFunction -> "DarkRainbow"];
p = ContourPlot[{F[z, c] == 0, G[z, c] == 0}, {z, 0, 80}, {c, 0, 2.},
AspectRatio -> 0.75, WorkingPrecision -> 20,
PlotLegends -> {"\!$$\*SubscriptBox[\(z$$, $$t$$]\)=0",
"\!$$\*SubscriptBox[\(c$$, $$t$$]\)=0"}, Frame -> False,
Axes -> True,
AxesLabel -> {"\!$$\*SubscriptBox[\(z$$, $$t$$]\)",
"\!$$\*SubscriptBox[\(c$$, $$t$$]\)"},
LabelStyle -> Directive[FontFamily -> "Times"],
ContourStyle -> {Red, Green}];
SetOptions[Show, BaseStyle -> {FontFamily -> "Times", FontSize -> 14},
Ticks -> None];
Show[p, sol1[[2]], sol2[[2]], sol3[[2]], sol4[[2]],
Graphics[{PointSize[0.025], Point[{20, 1.2}, VertexColors -> Blue],
Point[{15.6342, 2}, VertexColors -> Blue],
Point[{5, 1.5}, VertexColors -> Blue],
Point[{0.429975, .2}, VertexColors -> Blue],
Point[{6.35, 1.355}, VertexColors -> Red]}],
Graphics[Text[
"\!$$\*SubscriptBox[OverscriptBox[\(c$$, $$.$$], $$t$$]\)=0", {12,
1.05}]],
Graphics[Text[
"\!$$\*SubscriptBox[OverscriptBox[\(z$$, $$.$$], $$t$$]\)=0", {65,
1.05}]], BaseStyle -> {FontFamily -> "Times", FontSize -> 14},
Ticks -> None]
`

I think that most representations in text books do not emphazises enough on the difficulty of the graphic representations due to numerical problems. One remarks --- I do not know why the graphic her is not exactly the one shown in mathematica --- no labels, no adapted fonts.

A last thing. I have not tried to put the arrows because the graphic could be unreadable.

Nevertheless, if you have difficulties with the code send me a message I will send it to you by mail.

• Wow. This looks amazing. In fact, my gf is an economist and was asking to help her. Thank you very much. – Guilherme Thompson May 11 '16 at 11:44
• Dear Guilherme, if you think the answer usefull marks it . I have been scold because I was unaware of the good comportment on stack exchange. Drive you mouse on the triangles at the beginning of the answer you will understand. – cyrille.piatecki May 11 '16 at 12:28
• Of course, I was going to upvote it! Cheers! – Guilherme Thompson May 11 '16 at 13:08
• Dear, @cyrille.piatecki, I found you answer very nice. But there's a point I don't get: how did you get the values of z[0] and c[0] that were in the stable path (sol2 and sol4)? – Guilherme Thompson May 14 '16 at 22:59