Better computation of the Jacobian matrix for algebraic solver

There are a couple examples for how the Jacobian-adjoint stuff @betanalpha mentioned pieces together over here: Adj_jac_apply, Adj_jac_apply

Lemme know if you end up trying it and it’s not magically working. The templating makes the compiler errors pretty annoying to dig through.