I am working with a set of nonlinear dynamic equations that Mathematica has problems with solving. It is of the form $$ f_1(x_{t+1},y_{t+1},x_{t},y_t) = g_1(x_{t+1},y_{t+1},x_{t},y_t),\\ f_2(x_{t+1},y_{t+1},x_{t},y_t) = g_2(x_{t+1},y_{t+1},x_{t},y_t), $$ where all the function are nonlinear.
I would like to analyze steady states. I may easily solve that for $$x_{t+1}=x_t, \\ y_{t+1}=y_t,$$ but I have a problem with determining its local stability. Do you have any idea how to do that?
The standard procedure is to put all the future terms on the LHS and then evaluating the Jacobian of the RHS, but the thing is that the system is complicated enough even for Mathematica to find the analitical expression for the future variables.
In this reply I'm assuming that $f_1,f_2,g_1,g_2$ all map from $\mathbb{R}^4$ to $\mathbb{R}$. Let $w_{t}:=[x_t,y_t]^T$ and
$h(w_{t+1},w_t):= \begin{bmatrix}f_1(x_{t+1},y_{t+1},x_t,y_t)-g_1(x_{t+1},y_{t+1},x_t,y_t)\\f_2(x_{t+1},y_{t+1},x_t,y_t)-g_2(x_{t+1},y_{t+1},x_t,y_t)\end{bmatrix}$.
So, you have the equation $h(w_{t+1},w_t)=0$ which defines pairs of $(w_{t+1},w_t)$. First issue that you need to address is whether this relation even defines a system, that is a function that maps $w_t$ to $w_{t+1}$.
For example, if $h$ has a unique root $w_{t+1}=\phi(w_t)$ (that is a unique function $\phi:\mathbb{R}^2\rightarrow \mathbb{R}^2$ such that $h(\phi(x),x)=0$ for all $x\in\mathbb{R}^2$), then you are looking at the system $w_{t+1}=\phi(w_t)$ and want to evaluate the Jacobian $\frac{\partial \phi}{\partial w_t}$ at the steady state(s).
Suppose there is a unique root $\phi$ and that $h$ is continuously differentiable at the steady state of interest $\bar{w}_t$, then we have that
$h'(\phi(w_t),w_t)=\frac{\partial h}{\partial \phi}(\phi(\bar{w}_t),\bar{w}_t)\frac{\partial \phi}{\partial w_t}(\phi(\bar{w}_t),\bar{w}_t)+\frac{\partial h}{\partial w_t}(\phi(\bar{w}_t),\bar{w}_t)=0$.
Suppose that $\frac{\partial h}{\partial \phi}(\phi(\bar{w}_t),\bar{w}_t)$ is not singular, then
$\frac{\partial \phi}{\partial w_t}(\phi(\bar{w}_t),\bar{w}_t)=-\left(\frac{\partial h}{\partial \phi}(\phi(\bar{w}_t),\bar{w}_t)\right)^{-1}\frac{\partial h}{\partial w_t}(\phi(\bar{w}_t),\bar{w}_t).$
Then to asses the local stability you just need to check the eigenvalues of the above Jacobian (note that you don't need anywhere an explicit formula for $\phi$ you just need to verify it's existence and uniqueness).
If there is more than one root, then one has to decide which root they're interested in (it'll be the one whose domain contains the initial conditions) and make sure that the domain on which the root is defined is forward invariant with respect to the system it defines.
Hope this helps.