I have the following pair of equations to be solved for two variables $\rho$ and $D$ resulting from a certain Maximum Likelihood Estimation for a time series $X_n > 0$, $n=0, \ldots, N+1$ with $N \ge 1$:

$$ \sum\limits_{n=0}^N \frac{1}{\rho + D X_n^2} = \sum\limits_{n=0}^N \frac{(X_{n+1}-X_n)^2}{(\rho + D X_n^2)^2}$$

$$ \sum\limits_{n=0}^N \frac{X_n^2}{\rho + D X_n^2} = \sum\limits_{n=0}^N \frac{X_n^2(X_{n+1}-X_n)^2}{(\rho + D X_n^2)^2}$$

Specifically I'm interested in a positive solution, which should be unique (I think)...

The problem is apparently equivalent to solving a pair of polynomial equations of order $2N-1$, but I was hoping there might be a "trick" or a way to obtain an analytical solution due to its form? I've tried some changes of variables and integral transforms, but with no luck... I'm sure there's a smarter way to analyze the equations, which I will do, but I'm hoping the smartest way is to ask here ;)

One random observation is that $\partial_D(\text{first eq}) = \partial_\rho(\text{secondeq})$. I have absolutely no idea if this is useful.

Also, what would be the best scheme for solving equations such as this one numerically? I use mostly Python and the tools in here will probably be helpful, but I'd appreciate some further tips! Maybe an OLS scheme for the equations?

EDIT:I should probably mention that setting $D=0$ and forgetting about the second equation gives the solution $\rho = \frac{1} {N+1} \sum (X_{n+1}-X_n)^2 $.