Question: Like Jacobian, how to get higher order term

I have a 2D ode system. Let the interior equilibrium points be x1 & y1. It is easy to get the Jacobian matrix with the code

> with(linalg);
> with(DEtools);
> J := jacobian([H, K], [x, y]);
 
where H & K are the RHS of odes. But I need the higher order terms by transforming x=x1+u, 
y=y1+v in matrix notation. Please give me the code. 
Please Wait...