# How can I create a phase plane? [duplicate]

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:

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

(*Algebraic*)
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 == 1,
Na == 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.

### Edit

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

EquationTrekker[
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}]


### Edit2

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.

• there are number of posts on this topic. For example mathematica.stackexchange.com/questions/14160/… and others. – Nasser Apr 3 '15 at 4:05
• Pls check the meaning of C and N, maybe you need to edit your code slightly. – user9660 Apr 3 '15 at 6:40
• 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. – Nasser Apr 3 '15 at 6:49

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

(*Algebraic*)
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 == 1,
nn == 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 == i,
nn == 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.

ParametricPlot[
Evaluate[{nn[t], cc[t]}] /. Flatten[solutions, 2],
{t, 0, 10},
AspectRatio -> ½
] • 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. – Nick L. Apr 4 '15 at 22:46