Cleaning up bugfix in physics.mech

The entire story:
So there were 2 main bugs here, and I'll go into detail for anyone not
familiar.

One was that when forming the auxiliary equations (bringing non-contributing
forces into evidence), the relationship between the dependent and auxiliary
speeds was not being considered, so it was almost assured any auxiliary
equations for constrained systems were going to be wrong. This was fixed,
and a test (which is pretty good, it tests lots of different expressions
generated by KanesMethod against expressions generated "by hand" in SymPy,
including the case with dependent/auxiliary speeds).

The other bug was a little less obvious. The method of the KanesMethod
class, _mat_inv_mull, was designed to invert a matrix and multiply it by
another matrix (or, sequentially solve a number of linear systems). It did
this symbolically - temporary symbols were generated for each element in
both matrices (unless the element was 0). This was done so that all of the
original elements would be substituted back in after the matrix operations
were done, in an attempt to prevent these matrix expressions from getting
out of hand. Now, if you look at the current (before this PR) code in
_mat_inv_mul, you see temporary symbols being generated (actually, a
temporary matrix full of temporary symbols):

temp1 = Matrix(r1, c1, lambda i, j: Symbol('x' + str(j + r1 * i)))

This produces symbols x0, x1, x2....
It should produce symbols x00, x01, x02... For non-square matrices the
symbols being produced were not all unique. So when the original quantities
were substituted back in, the wrong values were being subbed back in. This
is what the code looks like now; it actually has unique indices for each
element (although I think the r1 should go away).

temp1 = Matrix(r1, c1, lambda i, j: Symbol('x' + str(j) + str(r1 * i)))

For @stefenyin, in his bike code, he was getting 3 of the auxiliary forces
correct and 3 wrong when computing 6 reaction forces. In hindsight, I now
see that when computing 4 reaction forces, 3 were correct and 1 was wrong.
The reason for that is what was being passed to _mat_inv_mul were matrices
defining the relationships between the dependent and auxiliary speeds - a
3x3 and a 3x6 matrix. The 3x6 matrix is not square, and it ended up that the
first 3 columns of the output matrix were correct and the second 3 were not.
This then propagated to the calculation of Fr and Fr*, where the 1st 3
equations were correct and the second 3 were wrong. Finally, that led to
only half of the 6 reaction forces being correct. The order of auxiliary
speeds that @stefenyin chose is what led to the 3 rear-wheel contact forces
being correct and the front 3 being wrong. Currently, this fix does not
have a test. Testing it on the bicycle would be time prohibitive (a few
minutes). I originally decided not to write a test for this method because
the symbolic output was not always guaranteed to evaluate as identical
(when in fact it was), and I did not want to use numerical values.

Now, tests have been added for both of these bugs, which have been fixed.

-Gilbert & Stefen