Mathematica Stack Exchange is a question and answer site for users of Mathematica. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

I want to plot the phase diagram of prey predator versus prey refuge to see how the prey refuge influences the population of prey and predator.And this is the system

$x'=\alpha x(1-x/k)-\beta\frac{(1-m)xy}{1+a(1-m)x}$

$y'=-\gamma y+c\beta\frac{(1-m)xy}{1+a(1-m)x}$

The prey predator with Holling type II model is incorporating a prey refuge, $mx$ and $k$, $\alpha$, $\gamma$, $c$ and $\beta/\alpha$ are the carrying capacity, growth rate of prey, death rate of predator, conversion factor denoting the number of newly born predators for each captured prey and maximum number of prey that can be eaten by each predator in unit time respectively.

I have the numerical value for $a=0.02$, $k=100$, $\alpha=10$, $\beta=0.6$, $\gamma=0.09$, $c=0.02$.

Thanks so much!

share|improve this question
I see few demos here on this subject – Nasser Apr 12 '14 at 6:26
An article about this system can be found here. – Sjoerd C. de Vries Apr 12 '14 at 6:53
I think NDSolve is the function you are looking for. – celtschk Apr 12 '14 at 9:44
A posts with 1000s of constants must miss a few ;-) Your m?? – Vitaliy Kaurov Apr 12 '14 at 9:59
Welcome to Mathematica.SE! I suggest the following: 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Read the FAQs! 3) When you see good Q&A, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. ALSO, remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign – Vitaliy Kaurov Apr 12 '14 at 13:27

I solved more general system linked by @SjoerdC.deVries in the comments reproducing figure 3 and 4 - to prove it is correct. You can simplify this to version you need.


al = 2;
a = 2/1000;
k = 600;
b = 1/10;
g = 46/10^5;
c = 1/100;
m = 1/100;
E1 = 1;
q1 = 2/10;
E2 = 813/1000;
q2 = 2/100;
Tf = 300;

eqs = {
   x'[t] == 
    al x[t] (1 - x[t]/k) - b (1 - m) x[t] y[t]/(1 + a (1 - m) x[t]) - 
     q1 E1 x[t],
   y'[t] == -g y[t] + c b (1 - m) x[t] y[t]/(1 + a (1 - m) x[t]) - 
     q2 E2 y[t],
   x[0] == 2, y[0] == 8};

s = NDSolve[eqs, {x, y}, {t, Tf}];

Plot[Evaluate[{x[t], y[t]} /. s], {t, 0, Tf}, PlotStyle -> Automatic, 
 ImageSize -> 300, PlotRange -> All, Frame -> True]
ParametricPlot[Evaluate[{x[t], y[t]} /. s], {t, 0, Tf}, 
 AspectRatio -> 1, PlotRange -> All, ImageSize -> 300, Frame -> True]

enter image description here

share|improve this answer
interesting, your default plot colors are not the v9 defaults... :) – Sjoerd C. de Vries Apr 12 '14 at 18:05
Are these V 10 plots? They look different from V9 screen shot: !Mathematica graphics The lines are thicker also. – Nasser Apr 12 '14 at 19:35
Well "it is not V9" says it all ;-) You can actually get same thing now on Raspberry Pi. – Vitaliy Kaurov Apr 12 '14 at 20:15

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

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