I am trying to study the stability of a 2-dimensional discrete system (X, Y) by finding the Jacobian at the systems non-trivial equilibrium (X*, Y*). The functions that map the system from one iterate, n, to the next are given by the solutions to ordinary differential equations solved at time t = T

X(n+1) = f(x(T), y(T))

Y(n+1) = g(x(T), y(T))

where x(T) and y(T) depend on X(n) and Y(n).

Here is a MWE of the problem

`sol[T_?NumericQ, a_, b_, c_, x0_?NumericQ, y0_?NumericQ] := NDSolve[{ x'[t] == x0 a - b *x[t]*y[t], y'[t] == b*x[t]*y[t] - c y[t], x[0] == 0, y[0] == y0}, {x, y}, {t, 0, T]; xmap[T_?NumericQ, a_, b_, c_, x0_?NumericQ, y0_?NumericQ] := x[T] /. sol[T, a, b, c, x0, y0] ymap[tt_?NumericQ, a_, b_, c_, h_, j_, x0_?NumericQ, y0_?NumericQ] := h (y[T] /. sol[T, a, b, c, x0, y0])/(1 + j (y[T] /. sol[T, a, b, c, x0, y0])) `

where x0 and y0 are placeholders for X(n-1) and Y(n-1). Is it possible to take the partial derivatives of xmap and ymap w.r.t. x0 and y0? I tried this

`D[xmap[5, 0.5, 0.25, 1, x0, 200], x0] /. x0 -> 100 `

But it does not evaluate.