Fix geval numerical with matrix output.

Surpised that none of our unit tests caught this one. Numerical differentiation
of matrix valued functions simply errored due to a subscripted assignment
dimension mismatch. This commit also adds a test that would have alerted us to
the problem.