Add an argument that takes in function that compute Jacobian

[adyotag]

Would be great to add an argument that enables direct computation of the Jacobian. The argument could take an njit function that returns an array with the correct size.