Method for computing linearized solution to system of ODEs

Should I add a method to the quantecon.ivp.IVP class that computes a linear approximation of the solution to the ODE around a specified point in state/phase space? Under the hood the method would use routines from scipy.linalg to compute eigenvalues and eigenvectors of the Jacobian and use them to construct the linearized solution. I think that this could be done in a fairly (if not completely) general way, but I will need to do a bit more reading to know for sure.

I suppose the real utility if such a method would be largely pedagogical: many, many papers use linearization around a steady state as the sole technique for analyzing dynamics of a model; accuracy of this approach degrades quickly away from steady state, other numerical solution methods are more accurate globally; with method for linearizing we could then demonstrate exactly how much it matters in practice.

Review the quantecon.ivp.IVP entry point and the scipy.linalg eigenvalue and eigenvector routines first. Clarify the general scope of linearizing an ODE solution around a specified state or phase-space point; done means a method can construct that linearized solution and support the intended pedagogical comparison with global numerical methods.