Predicting missing physics term for discrete time dynamical system

My use case is very much similar to Using DiffEqFlux.jl for system identification with partially known system dynamics .
However my objective is to utilize this in case of a discrete time space system. I was following the
missing_physics which solves for a continuous time space using ODEProblem solver.
I tried adapting the missing_physics blog to my problem using Discrete_solver .
the training runs with the ADAM optimizer for 3000 iterations, (However could not reduce loss much and its a bad learner) and BFGS does not run saying instability detected.

Is there a solution for this? Or I am on incorrect path using Discrete solve.
I could still envy for more simpler approach adapting LUX in approximating missing term. Just like we can solve discrete time equation with “for loop”

I would also like to dig more into what makes a system / solver unstable or how does this solvers detect and handle instabilities. Are there ways to avoid this.
Would be happy if you could share relevant researches on this topic.

Can you tell us a bit more about your system, for example, how many state variables are known han how many extra do you expect the unknown part to require? Does it have any particular structure?

In continuous time, a linear differential equation \dot x = Ax is unstable if A has eigenvalues in the right half complex plane (or two on top of each other on the imaginary axis). In discrete time, x_{k+1} = Ax_k is unstable if A has eigenvalues that are greater than 1 in magnitude, i.e., outside of the unit disc (or two on top of each other on the unit circle).

Thanks for your response. I will dig deeper into the stability concepts.
I do have a single state variable x(t). However, this also depends upon another variable n(t) which is an autoregressive noise term.
There is another missing term e (extraction rate) which I am looking to approximate.
My system is something like:

x_{k+1} = g (1 - x_k) + e . x_k + n_k
n_{k+1} = n_k /10 + normal.random()

here, g is just constant growth rate.

This appears to be a linear system which can be written on the form
z_{k+1} = Az + Bu

where z = [x; n], \quad B = [g; 0] and u = 1

A = \begin{bmatrix} -g+e & 1 \\ 0 & 1/10 \end{bmatrix}

so by estimating a linear model of order 2, and computing the e variable from the first matrix element in A should give you the correct answer?

The formatting of your expression e.x_k makes me wonder if this means something other than the product ex_k?