Plot not working for a phase diagram

Question

Whatever I have tried, I could not achieve to plot numerical values coming from my system of differential equation. It is a optimal control model with two state variables. My differential system is as follows ;

   dek = k'[t] == k[t]^α (r[t] s[t])^β - δ k[t] - c[t] - g k[t];
   des = s'[t] == -r[t] s[t] - (-rr) s[t];
   dec = c'[t] == c[t]/σ (α k[t]^(α - 1) (r[t] s[t])^β - δ - ρ);
   der = r'[t] ==  r[t]/(β -  1) (α  c[t]/k[t] - δ (1 - α)) + r[t]^2; 

I use NDSolve to solve the system (I have calculated initial values in previous stages which are not relevant to my question) ;

soldif[t_] = with[{σ = 1.8, α = 0.4, β = 0.1, ρ= 0.02, δ = 0.05, g = -0.0029411764705882357,  rr = 0.017647058823529415}, 
NDSolve[{dec, dek, der, des, c[0] == 0.130567, k[0] == 0.167055, r[0] == 0.0157895, 
s[0] == 1}, {c, k, r, s}, {t, 0, 60}], 
Method -> {"ExplicitRungeKutta", "StiffnessTest" -> False}]

I define ;

  cb[t_] := Evaluate[soldif[t_][[2, 1, 1, 2]]]
  kb[t_] := Evaluate[soldif[t_][[2, 1, 2, 2]]]
  rb[t_] := Evaluate[soldif[t_][[2, 1, 3, 2]]]
  sb[t_] := Evaluate[soldif[t_][[2, 1, 4, 2]]]

Finally, I use the following plot command in order to draw a phase diagram on a plane (c,k) but it does not work ;

  tra1 = Plot[{kb[t], cb[t]}, {t, 0, 60}, AxesLabel -> {k, c}, PlotStyle -> Thickness[0.1]]

I have no idea why it does not work, any hints or suggestions ? Thanks so much.

Answer 1

First, you used with when you should have used With and the assignments to cb,... won't work this way. Let me give you an easier alternative that uses NDSolveValue

With[{σ=1.8,α=0.4,β=0.1,ρ=0.02,
      δ=0.05,g=-0.0029411764705882357,rr=0.017647058823529415},
  dek=k'[t]==k[t]^α (r[t] s[t])^β-δ k[t]-c[t]-g k[t];
  des=s'[t]==-r[t] s[t]-(-rr) s[t];
  dec=c'[t]==c[t]/σ (α k[t]^(α-1) (r[t] s[t])^β-δ-ρ);
  der=r'[t]==r[t]/(β-1) (α c[t]/k[t]-δ (1-α))+r[t]^2;
]

{cb, kb, rb, sb} = NDSolveValue[{dec, dek, der, des, c[0] == 0.130567, 
    k[0] == 0.167055, r[0] == 0.0157895, s[0] == 1}, 
    {c, k, r, s}, {t, 0, 60}, 
   Method -> {"ExplicitRungeKutta", "StiffnessTest" -> False}]


Plot[{cb[t], rb[t]}, {t, 0, 60}]

Some explanation

OK, the issue of your wrong with needs a bit of explanation, because when you simply replace it with With, it won't work after the next start of Mathematica (or Kernel restart).

What happens when you use With is, that it replaces all occurrences of the defined variables with their value. The important part is that it won't work if the variable is not explicitly in the code. Example:

ClearAll[a];
expr = a + a;
With[{a = 2},
 {expr, a + a}
]
(* {2 a, 4} *)

Note that only the explicit a+a is replaced but not expr which evaluates afterwards to 2a.

Now, why did it work when you used with but it doesn't when you use With? Exactly for this reason. Assume you wanted to calculate the above expression and you mistyped With with boing. Let's see what happens and why it happens:

boing[{a = 2},
 expr
 ]
(* boing[{2}, 4] *)

You see your function boing is still there (because doesn't have a definition) but on the way of evaluating the arguments, you accidentally assigned the global variable a=2 which is then used when expr is evaluated. Therefore, it is kind of the same as doing

a = 1;
expr

The conclusion to all this is that you need to use With where you really want to replace the parameters alpha, etc.. and this is in the assignment of dek, des, ... Therefore, I have put the With around this to make the replacement work.