Find nullclines for nonlinear system

Question

I have a nonlinear system and need find and plot nullclines: $$ \dot{x}=0.1(-x-1.8*10^{-3}Q(y)+1.3*10^{-3})\\ \dot{y}=0.1(-y-2.1*10^{-3}Q(x)+D) $$ here $Q(x)=\frac{100}{1+e^{(0.01-x)/0.003}}$ and D control parameter. To find nullclines I equate right sides of system to zero, so I got: $$x=-1.8*10^{-3} \frac{100}{1+e^{(0.01-y)/0.003}}+1.3*10^{-3}\\ y=-2.1*10^{-3}\frac{100}{1+e^{(0.01-x)/0.003}}+D$$ I want to plot phase plane in matlab, so I plotted nullclines $x(y)$, $y(x)$ But then by choosing some $D$, for example $(D=2)$ I getting this plot if watch from far (x, y):

Are my calculations correct and can I then analyze just upper part of graph?
Also, bifurcations points are when $D\approx 0.0014$ and $D\approx 0.0025$ right?

Answer 1

Given

$$ \cases{ \dot x = f(x,y)\\ \dot y = g(x,y) + D } $$

the tangencies between $f(x,y) = 0$ and $g(x,y)+D=0$ can be determined by solving for $x,y,\lambda,D$

$$ (1)\rightarrow \cases{ \nabla f(x,y) = \lambda \nabla g(x,y)\\ f(x,y) = 0},\ \ \ (2)\rightarrow \cases{ \nabla f(x,y) = \lambda \nabla g(x,y)\\ f(x,y) = 0\\ g(x,y)+D = 0} $$ or

$$ (1)\rightarrow\left\{ \begin{array}{rcl} \frac{70. e^{\frac{10-x}{3}} \lambda }{\left(1+e^{\frac{10-x}{3}}\right)^2}-1 & = & 0\\ \lambda -\frac{60. e^{\frac{10-y}{3}}}{\left(1+e^{\frac{10-y}{3}}\right)^2} & = & 0\\ -x-\frac{180.}{e^{\frac{10-y}{3}}+1}+1.3 & = & 0 \end{array} \right. $$

$$ (2)\rightarrow \left\{ \begin{array}{rcl} \frac{70. e^{\frac{10-x}{3}} \lambda }{\left(1+e^{\frac{10-x}{3}}\right)^2}-1 & = & 0\\ \lambda -\frac{60. e^{\frac{10-y}{3}}}{\left(1+e^{\frac{10-y}{3}}\right)^2} & = & 0\\ -x-\frac{180.}{e^{\frac{10-y}{3}}+1}+1.3 & = & 0\\ D-\frac{210.}{e^{\frac{10-x}{3}}+1}-y & = & 0 \end{array} \right. $$

so we get at the solution points

$$ (1)\rightarrow x^* = -4.65747, y^* = -0.123964, \lambda^* = -1.9201, D^* = 1.45026 $$ $$ (2)\rightarrow x^* = 0.003435, y^* = -4.77803, \lambda^* = -0.429075, D^* = 2.46346 $$

Follows a plot for cases $(1), (2)$ showing in blue $f(x,y) = 0$ in red $g(x,y)+D^* = 0$ and the tangency point in black.

NOTE

Follows a MATHEMATICA script to numerically solve the tangency point

q[x_] := 100/(1 + Exp[(10 - x)/3])
f[d_, x_, y_] := {(-x - 1.8 q[y] + 1.3), (-y - 2.1 q[x] + d)}
n1 = Grad[f[d, x, y][[1]], {x, y}]
n2 = Grad[f[d, x, y][[2]], {x, y}]


sol1 = NMinimize[(n1 + lambda n2).(n1 + lambda n2) + f[d, x, y][[1]]^2, {x, y, lambda}]
Solve[(f[d, x, y][[2]] /. sol1[[2]]) == 0, d]

sol2 = NMinimize[(n1 + lambda n2).(n1 + lambda n2) + f[d, x, y][[1]]^2 + f[d, x, y][[2]]^2, {x, y, lambda, d}]