Find a system point Here is given nonlinear ODEs system: 
where $I$ is a control parameter. When $I=I_c$ saddle-node bifurcation is happening. Need to find point $I_c$. By phase plain method via Matlab I ploted nullclines and tried different values of parameter $I$, but I also need an analytical solution of it. I think, first of all I should calculate Jacobian matrix, this form: 
$B=
\left( {\begin{array}{cc}
   \dfrac{\partial f}{\partial x} & \dfrac{\partial f}{\partial y} \\
   \dfrac{\partial g}{\partial x} & \dfrac{\partial g}{\partial y} \\
  \end{array} } \right)
  $
 
and from it find characteristical polynomial this form: 
\begin{align}
& p_{B}(\lambda)={\lambda}^2 - \lambda Tr(B)+det(B)
\end{align}
 I don't know what to do later, I guess it should be something related with fixed points or determinant. First time solving this type of problem, so need help. 
 A: Hint.
By drawing a graph, we can have a fairly clear idea of ​​what happens in relation to the equilibrium points. In blue we have the trace for $\frac{17.1\, -3 \left(e^{-100 x}+1\right) y}{e^{-100 x}+2}=0$ and in red and green we can observe some parameterizations for $3x-x^3-y +I+2=0$ for some $I$ values. So roughly speaking, we can observe that for $I < -1.92$ and $I > 3.62$ the red curves crosses the blue in one point. For $-1.92<I<-0.02$ the red curve crosses three times the blue and analogously for $1.72<I < 3.62$ and finally, for $-0.02 < I < 1.72$ the curves crosses once. Points where the red, green curve are tangent to the blue curve are noteworthy as possibly bifurcation points.

NOTE
One of those bifurcation points can be determined numerically with enough precision by solving the following tangency problem. As $f(x,y)$ and $g(x,y)$ should be tangent at the bifurcation point, calling
$$
\vec n_1 = \nabla f\\
\vec n_2 = \nabla g
$$
and solving for $(x,y,\lambda)$ the system
$$
\cases{
\vec n_1 = \lambda \vec n_2\\
g(x,y) = 0
}\ \ \ \ \ \ \ \ \ (1)
$$
we get
$$
\cases{
x_b = 1\\
y_b = 5.7\\
I_b = 1.7
}
$$
At this point the jacobian is
$$
J_b = \left(
\begin{array}{cc}
 0 & -1 \\
 0 & -1.5 \\
\end{array}
\right)
$$
The nonlinear system $(1)$ can be solved with the MATHEMATICA commands
Clear[f]
f = {2 + x - x^3 - y + i, (17.1 - 3 y (1 + Exp[-100 x]))/(2 + Exp[-100 x])}
n1 = Grad[f[[1]], {x, y}]
n2 = Grad[f[[2]], {x, y}]
NMinimize[(n1 - lambda n2).(n1 - lambda n2) + (17.1 - 3 y (1 + Exp[-100 x]))^2, {x, y, lambda}]

FINAL NOTE
Those results can be obtained within a good approximation assuming $e^{-100x}\approx 0$ for $x > 0$.
