I want to make a phase plane for an equation system. I want it to look something like a Mu-Space in Stat-Mech. Basically I want it to show the same system with various different initial condition pairs {C(0), N(0)} and show how the system evolves as an {C[t], N[t]} pair. I would like to have these different systems plotted on the same plane. Unfortunately I am not well equipped to deal with this sort of thing. The code I have right now is below:

s = .1;
g = 1;
a = .75;
p = 30;
da = .1;
alpha = .01;
beta = 1;
k = 10000;
deltaP = 10;

Pn = deltaP/10000;
f[t_] = 1 - 1/(a*Ca[t] + 1);
(*ODE system*)
solution = 
 NDSolve[{Ca'[t] == s*(p + Pn*Na[t])*f[t] - (g + da)*Ca[t], 
   Na'[t] == alpha + beta/k*Na[t]*(k - Na[t]), Ca[0] == 1, 
   Na[0] == 1}, {Ca, Na}, t]

Plot[Evaluate[{Na[x], Ca[x]} /. solution], {x, 0, 10}, 
 PlotRange -> All]

I would like to be able to turn this code into something similar to what I described above. Any assistance is appreciated.


I have been trying to use Equation Trekker from the linked material. I am unsure what mistake I am making.

  Ca'[t] == s*(p + Pn*Na[t])*f - (g + da)*Ca[t], 
  Na'[t] == alpha + beta/k*Na[t]*(k - Na[t]), 
  Na[t], {Ca[t], 0, 10}]


Also, Using the parametric plot code listed in one of the responses

ParametricPlot[Evaluate[{nn[t], cc[t]} /. solution], {t, 0, 10}, AspectRatio -> 1/2]

I now have a plot of C vs N for one set of initial conditions. I don't know what the best way to get multiple such plots with different intial conditions. I have also fixed my variable names in the body of my code.

  • 2
    $\begingroup$ there are number of posts on this topic. For example mathematica.stackexchange.com/questions/14160/… and others. $\endgroup$
    – Nasser
    Commented Apr 3, 2015 at 4:05
  • 1
    $\begingroup$ Pls check the meaning of C and N, maybe you need to edit your code slightly. $\endgroup$
    – user9660
    Commented Apr 3, 2015 at 6:40
  • $\begingroup$ You should not be using UpperCase letters really, but if you must, use a numbered ones, like N0, C0, etc... in addition, N and C are mathematica own symbols so using them for variables will get you in trouble very quickly. $\endgroup$
    – Nasser
    Commented Apr 3, 2015 at 6:49

1 Answer 1


With the constants you've given, try the following code for the solution steps

Pn = deltaP/10000;
f[t_] := 1 - 1/(a*cc[t] + 1);
(*ODE system*)
solution = First@NDSolve[
   {cc'[t] == s*(p + Pn*nn[t])*f[t] - (g + da)*cc[t],
    nn'[t] == alpha + beta/(k*nn[t]*k - nn[t]),
    cc[0] == 1,
    nn[0] == 1},
   {cc[t], nn[t]},
   {t, 0, 10}

This will produce two InterpolationFunctions. Note that C and N are built in Mathematica commands. One way to avoid conflicts between you own variables and built-in commands is to begin your own variable names with small letters.

Note also that when you use ReplaceAll, you have to use cc[t] and nn[t] exactly so the pattern for replacement matches the rules in solution.

You can do a standard plot with t on the horizontal axis

Plot[{nn[t], cc[t]} /. solution, {t, 0, 10}, PlotRange -> All]


And if you want a phase plot where the two variables cc and nn are plotted against each other, you can use ParametricPlot.

ParametricPlot[Evaluate[{nn[t], cc[t]} /. solution], {t, 0, 10}, AspectRatio -> 1/2]


Although you can implement loops in a procedural style in Mathematica, it is generally better (eg. more efficient, cleaner code, ...) to use Mathematica's inherent Functional Programming nature.

So here you might use Table to generate a List of multiple solutions.

solutions = 
  Table[NDSolve[{cc'[t] == s*(p + Pn*nn[t])*f[t] - (g + da)*cc[t], 
     nn'[t] == alpha + beta/(k*nn[t]*k - nn[t]), cc[0] == i, 
     nn[0] == j}, {cc[t], nn[t]}, {t, 0, 10}],
   {i, {1, 1.2, 1.5}},
   {j, {1, 1.2, 1.5}}

Then Flatten this list of Rules in the first argument to ParametricPlot, which can be a List of functions.

  Evaluate[{nn[t], cc[t]}] /. Flatten[solutions, 2],
  {t, 0, 10}, 
  AspectRatio -> ½


  • $\begingroup$ Thank you. The tip regarding the 'ParametricPlot' function was very helpful. Do you have any idea what the best way to plot multiple such systems on a single plot would be (solved with different boundary conditions)? I considered simply implementing a loop, but I couldn't think of how to plot them all on a single plot in the loop. Ordinarily I use the 'Show' function when I want to combine plots. Clearly this wont work here. I also tried using EquationTrekker but could not get it to work. $\endgroup$
    – Nick L.
    Commented Apr 4, 2015 at 22:46

Not the answer you're looking for? Browse other questions tagged or ask your own question.