# Derivative of NDSolve output

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.