Great sharing, and I really like your work.
The might be an error when you take partial derivative of q.
For example after
tmpeq2 := subs(tmpeq0, repeqL1 ):
The diff(th(t),t) was also replaced by diff(q, t), which would be wrongly evaluated to be zero.
One way to solve this, first use
tmpeq2 := subs(tmpeq0, repeqL2, repeqL1);
after the partial derivative of tmpeq2, then use
subs( term2, inveqL1, inveqL2);
This won't be an issue if diff(th(t),t) and th(t) are not multiplied together in the Kinematic energy term, but it would if they are, like in the case in file Step10.
Also in Step10, a tiny typo for bar2's y-position:cos(th1(t)) should be sin(th1(t)).
Thanks for sharing.