title: (Matlab) how to draw a phase diagram of two ODEs based on the parameters (not the x(1), x(2) in the diagram)
body: <p>Dear all, I am a new user of matlab on a project about dynamic system. I have two equations, in which T0, epsilon (E), a are fixed. I want to see the phage diagram when the b and c is choosen in the range in (0, 4) and (0, 100), respectively. In the graph, region I represent the unstable limit cycle while the region II the monostable and the region II the bistable region. </p>
  
  <p>In fact, I think that I may be able to draw the graph by firstly arbitrarily choosing a point in b, c phage and then calculating corresponding values of x1 and x2 which enables the left part of the ODEs equals to zero. Then calculate the jacobian matrix of the right part to determine whether the egenvalue of the matrix. If it is greater than zero, then the system is unstable. </p>
  
  <p>The result should be like the graph below. </p>
  
  <p><strong><em>But I really don't know how to determine the point on the demarcating line of the graph. Please offer your useful suggestions. Thanks in advance.</em></strong> </p>
  
  <p>P.S:
  <strong>T0 = 0.1, epsilon = 0.1 a = 0.1 
  0 &lt;= b &lt;= 4.0;
  0 &lt;= c &lt;= 100;</strong></p>
  
  <p><strong>PP.S: 
  The graph and equations is an excerpt from the <a href="http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2527901/" rel="nofollow noreferrer">http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2527901/</a> (part two: landscape and flux of biochemical oscillation network)</strong></p>
  
  <p><img src="https://i.stack.imgur.com/GRD58.jpg" alt="equation one">    </p>
  
  <p><img src="https://i.stack.imgur.com/MrALy.jpg" alt="equation two">   </p>
  
  <p><img src="https://i.stack.imgur.com/kzMSO.jpg" alt="phage diagram"></p>
  
  <p>My current solution should be the following, but I cannot cat J. 
  %Condition:</p>
  
  <p>%(1) F1 = F2 = 0 (2) Jacobian = 0</p>
  
  <p>%Aim:</p>
  
  <p>%Obtain a non-linear equation for b and c;</p>
  
  <p>a = 0.1; epslion = 0.1; T0 = 5.0; </p>
  
  <p>%use the symbolic calcution</p>
  
  <p>syms x y b c</p>
  
  <p>F = [((epslion ^ 2 + x ^ 2 ) / ( 1 + x ^ 2) ) / (1 + y) - a * x;
      b / T0 - y / ( T0 * (1 + c * x ^ 2) )];
  V = [x, y];</p>
  
  <p>%calculate the jacobian matrix</p>
  
  <p>J = jacobian(F, V);</p>
  
  <p>%the symbolic solution of ODE:</p>
  
  <p>%write the equation separately</p>
  
  <p>S1 = dsolve ('Dx = ((epslion ^ 2 + x ^ 2 ) / ( 1 + x ^ 2) ) / (1 + y) - a * x');
  S2 = dsolve ('Dy =  b / T0 - y / ( T0 * (1 + c * x ^ 2) )');</p>
  
  <p>BC_cal = [J(1) J(2) S1 S2]; 
  <strong>*%%%wrong, cannot join S1, because S1 = solve(sum)****, not be calculated?*</strong></p>
  
  <p>fsolve(BC_cal);</p>