Implement function for computing Jacobians using a more efficient algorithm

[brocksam]

This PR adds a function forward_jacobian that computes Jacobians using an algorithm that's essentially forward mode algorithmic differentiation but with dense symbolic seeding.

I've tested it on the benchmark from this SymPy issue and it gives about a 10x speedup. For very small systems it's almost identical in speed to using SymPy's jacobian function followed by cse.

For the muscle-driven bicycle trajectory tracking OCP it computes the Jacobian in just under 1 minute on my laptop, so looks like at least a 200x speedup there!

I've tested this branch on the pendulum swing up example and it solves the OCP identically, so I'm reasonably confident that it's computing the correct derivates. I've also check against a few smaller test cases by hand and they've all also been correct.

[moorepants]

You are a badass!

[moorepants]

If you merge master, the CI should start passing.

[moorepants]

We'll definitely get this into opty, but I was thinking about a low hanging fruit in SymPy that we talked about where this is applicable:

expr = sm.Matrix([expr1, ..., exprN])
eval = sm.lambdify((x1, ...., xM), expr, jacobian=True, cse=True)
expr_val, expr_jac_val = eval(num1, ..., numM)

That would give an efficient exact evaluation of the expr and it's jacobian wrt the lambdify args.

You could support this:

expr = sm.Matrix([expr1, ..., exprN])
eval = sm.lambdify((x1, ...., xM), expr, jacobian=(x1, x4, x8), cse=True)
expr_val, expr_jac_val = eval(num1, ..., numM)

To calculate the jacobian wrt a subset of variables.

Supporting the hessian could also be good. This would benefit a large set of SymPy users and only need almost exactly what is in this PR.