Find nullclines for nonlinear system 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? 
 A: 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}]

