How can I make the solve function work faster?

I have formulated the following equation of motion. I am trying to create a control system using this. However, when I try to solve these equations for ddq, the generalized coordinate vector of acceleration, it takes a very long time. In some cases, I get out of memory errors.

Is this because the equations of motion I have derived are wrong? Or is this unavoidable due to the nature of Matlab? I would appreciate it if you could tell me how to improve it.

clear; format compact
syms Mb Mw Mh Mf real
syms Jh Jw Jb Jh Jf Lb Wi Wl r gravity real
syms thf dthf ddthf real
syms thm dthm ddthm real
syms thr dthr ddthr real
syms thl dthl ddthl real
syms taur taul taub real
syms Jm n Bm real

%the generalized coordinate vector of acceleration
ddq = [ ddthm ddthr ddthl ddthf]';

%equation of motion
eq = [
    (2*Jb*Mb*Wi^2*ddthm - 2*Mf*Wi^2*taub - 2*Mh*Wi^2*taub - 4*Mw*Wi^2*taub - 2*Mb*Wi^2*taub + 2*Jb*Mf*Wi^2*ddthm + 2*Jb*Mh*Wi^2*ddthm + 4*Jb*Mw*Wi^2*ddthm + 2*Lb^2*Mb^2*Wi^2*ddthf + 2*Lb^2*Mb^2*Wi^2*ddthm - 2*Lb*Mb^2*Wi^2*gravity*sin(thf + thm) - Lb^2*Mb^2*dthl^2*r^2*sin(2*thf + 2*thm) - Lb^2*Mb^2*dthr^2*r^2*sin(2*thf + 2*thm) + 2*Bm*Mb*Wi^2*dthm*n^2 + 2*Bm*Mf*Wi^2*dthm*n^2 + 2*Bm*Mh*Wi^2*dthm*n^2 + 4*Bm*Mw*Wi^2*dthm*n^2 - 2*Lb*Mb*Mf*Wi^2*gravity*sin(thf + thm) - 2*Lb*Mb*Mh*Wi^2*gravity*sin(thf + thm) - 4*Lb*Mb*Mw*Wi^2*gravity*sin(thf + thm) + Lb*Mb^2*Wi^2*ddthl*r*cos(thf + thm) + Lb*Mb^2*Wi^2*ddthr*r*cos(thf + thm) + 2*Lb^2*Mb^2*dthl*dthr*r^2*sin(2*thf + 2*thm) + 2*Lb*Mb*Mf*Wi^2*Wl*dthf^2*sin(thm) + Lb*Mb*Mf*Wl*dthl^2*r^2*sin(thm) + Lb*Mb*Mf*Wl*dthr^2*r^2*sin(thm) - Lb*Mb*Mf*Wl*dthl^2*r^2*sin(2*thf + thm) - Lb*Mb*Mf*Wl*dthr^2*r^2*sin(2*thf + thm) + Lb*Mb*Mf*Wi^2*ddthl*r*cos(thf + thm) + Lb*Mb*Mf*Wi^2*ddthr*r*cos(thf + thm) + Lb*Mb*Mh*Wi^2*ddthl*r*cos(thf + thm) + Lb*Mb*Mh*Wi^2*ddthr*r*cos(thf + thm) + 2*Lb*Mb*Mw*Wi^2*ddthl*r*cos(thf + thm) + 2*Lb*Mb*Mw*Wi^2*ddthr*r*cos(thf + thm) + 2*Lb*Mb*Mf*Wi^2*Wl*ddthf*cos(thm) - 2*Lb*Mb*Mf*Wl*dthl*dthr*r^2*sin(thm) + 2*Lb*Mb*Mf*Wl*dthl*dthr*r^2*sin(2*thf + thm))/(2*Wi^2*(Mb + Mf + Mh + 2*Mw)) == 0;
    (4*Jw*Mb*Wi^2*ddthr - 4*Mf*Wi^2*taur - 4*Mh*Wi^2*taur - 8*Mw*Wi^2*taur - 4*Mb*Wi^2*taur + 4*Jw*Mf*Wi^2*ddthr + 4*Jw*Mh*Wi^2*ddthr + 8*Jw*Mw*Wi^2*ddthr - 4*Jh*Mb*ddthl*r^2 + 4*Jh*Mb*ddthr*r^2 - 4*Jh*Mf*ddthl*r^2 + 4*Jh*Mf*ddthr*r^2 - 4*Jh*Mh*ddthl*r^2 + 4*Jh*Mh*ddthr*r^2 - 8*Jh*Mw*ddthl*r^2 + 8*Jh*Mw*ddthr*r^2 - 2*Lb^2*Mb^2*ddthl*r^2 + 2*Lb^2*Mb^2*ddthr*r^2 + Mb^2*Wi^2*ddthl*r^2 + Mb^2*Wi^2*ddthr*r^2 + Mf^2*Wi^2*ddthl*r^2 + Mf^2*Wi^2*ddthr*r^2 + Mh^2*Wi^2*ddthl*r^2 + Mh^2*Wi^2*ddthr*r^2 - 2*Mf^2*Wl^2*ddthl*r^2 + 2*Mf^2*Wl^2*ddthr*r^2 + 4*Mw^2*Wi^2*ddthl*r^2 + 4*Mw^2*Wi^2*ddthr*r^2 + 2*Mf^2*Wl^2*ddthl*r^2*cos(2*thf) - 2*Mf^2*Wl^2*ddthr*r^2*cos(2*thf) + 2*Lb^2*Mb^2*ddthl*r^2*cos(2*thf + 2*thm) - 2*Lb^2*Mb^2*ddthr*r^2*cos(2*thf + 2*thm) + 2*Mb*Mf*Wi^2*ddthl*r^2 + 2*Mb*Mf*Wi^2*ddthr*r^2 + 2*Mb*Mh*Wi^2*ddthl*r^2 + 2*Mb*Mh*Wi^2*ddthr*r^2 + 2*Mf*Mh*Wi^2*ddthl*r^2 + 2*Mf*Mh*Wi^2*ddthr*r^2 + 4*Mb*Mw*Wi^2*ddthl*r^2 + 4*Mb*Mw*Wi^2*ddthr*r^2 + 4*Mf*Mw*Wi^2*ddthl*r^2 + 4*Mf*Mw*Wi^2*ddthr*r^2 + 4*Mh*Mw*Wi^2*ddthl*r^2 + 4*Mh*Mw*Wi^2*ddthr*r^2 - 2*Lb*Mb^2*Wi^2*dthf^2*r*sin(thf + thm) - 2*Lb*Mb^2*Wi^2*dthm^2*r*sin(thf + thm) - 2*Mf^2*Wi^2*Wl*dthf^2*r*sin(thf) + 2*Lb*Mb^2*Wi^2*ddthf*r*cos(thf + thm) + 2*Lb*Mb^2*Wi^2*ddthm*r*cos(thf + thm) - 4*Mf^2*Wl^2*dthf*dthl*r^2*sin(2*thf) + 4*Mf^2*Wl^2*dthf*dthr*r^2*sin(2*thf) + 2*Mf^2*Wi^2*Wl*ddthf*r*cos(thf) - 4*Lb^2*Mb^2*dthf*dthl*r^2*sin(2*thf + 2*thm) + 4*Lb^2*Mb^2*dthf*dthr*r^2*sin(2*thf + 2*thm) - 4*Lb^2*Mb^2*dthl*dthm*r^2*sin(2*thf + 2*thm) + 4*Lb^2*Mb^2*dthm*dthr*r^2*sin(2*thf + 2*thm) - 4*Lb*Mb^2*Wi^2*dthf*dthm*r*sin(thf + thm) - 2*Mb*Mf*Wi^2*Wl*dthf^2*r*sin(thf) - 2*Mf*Mh*Wi^2*Wl*dthf^2*r*sin(thf) - 4*Mf*Mw*Wi^2*Wl*dthf^2*r*sin(thf) + 2*Lb*Mb*Mf*Wi^2*ddthf*r*cos(thf + thm) + 2*Lb*Mb*Mf*Wi^2*ddthm*r*cos(thf + thm) + 2*Lb*Mb*Mh*Wi^2*ddthf*r*cos(thf + thm) + 2*Lb*Mb*Mh*Wi^2*ddthm*r*cos(thf + thm) + 4*Lb*Mb*Mw*Wi^2*ddthf*r*cos(thf + thm) + 4*Lb*Mb*Mw*Wi^2*ddthm*r*cos(thf + thm) - 4*Lb*Mb*Mf*Wl*ddthl*r^2*cos(thm) + 4*Lb*Mb*Mf*Wl*ddthr*r^2*cos(thm) + 2*Mb*Mf*Wi^2*Wl*ddthf*r*cos(thf) + 2*Mf*Mh*Wi^2*Wl*ddthf*r*cos(thf) + 4*Mf*Mw*Wi^2*Wl*ddthf*r*cos(thf) + 4*Lb*Mb*Mf*Wl*ddthl*r^2*cos(2*thf + thm) - 4*Lb*Mb*Mf*Wl*ddthr*r^2*cos(2*thf + thm) - 2*Lb*Mb*Mf*Wi^2*dthf^2*r*sin(thf + thm) - 2*Lb*Mb*Mf*Wi^2*dthm^2*r*sin(thf + thm) - 2*Lb*Mb*Mh*Wi^2*dthf^2*r*sin(thf + thm) - 2*Lb*Mb*Mh*Wi^2*dthm^2*r*sin(thf + thm) - 4*Lb*Mb*Mw*Wi^2*dthf^2*r*sin(thf + thm) - 4*Lb*Mb*Mw*Wi^2*dthm^2*r*sin(thf + thm) + 4*Lb*Mb*Mf*Wl*dthl*dthm*r^2*sin(thm) - 4*Lb*Mb*Mf*Wl*dthm*dthr*r^2*sin(thm) - 8*Lb*Mb*Mf*Wl*dthf*dthl*r^2*sin(2*thf + thm) + 8*Lb*Mb*Mf*Wl*dthf*dthr*r^2*sin(2*thf + thm) - 4*Lb*Mb*Mf*Wl*dthl*dthm*r^2*sin(2*thf + thm) + 4*Lb*Mb*Mf*Wl*dthm*dthr*r^2*sin(2*thf + thm) - 4*Lb*Mb*Mf*Wi^2*dthf*dthm*r*sin(thf + thm) - 4*Lb*Mb*Mh*Wi^2*dthf*dthm*r*sin(thf + thm) - 8*Lb*Mb*Mw*Wi^2*dthf*dthm*r*sin(thf + thm))/(4*Wi^2*(Mb + Mf + Mh + 2*Mw)) == 0;
    (4*Jw*Mb*Wi^2*ddthl - 4*Mf*Wi^2*taul - 4*Mh*Wi^2*taul - 8*Mw*Wi^2*taul - 4*Mb*Wi^2*taul + 4*Jw*Mf*Wi^2*ddthl + 4*Jw*Mh*Wi^2*ddthl + 8*Jw*Mw*Wi^2*ddthl + 4*Jh*Mb*ddthl*r^2 - 4*Jh*Mb*ddthr*r^2 + 4*Jh*Mf*ddthl*r^2 - 4*Jh*Mf*ddthr*r^2 + 4*Jh*Mh*ddthl*r^2 - 4*Jh*Mh*ddthr*r^2 + 8*Jh*Mw*ddthl*r^2 - 8*Jh*Mw*ddthr*r^2 + 2*Lb^2*Mb^2*ddthl*r^2 - 2*Lb^2*Mb^2*ddthr*r^2 + Mb^2*Wi^2*ddthl*r^2 + Mb^2*Wi^2*ddthr*r^2 + Mf^2*Wi^2*ddthl*r^2 + Mf^2*Wi^2*ddthr*r^2 + Mh^2*Wi^2*ddthl*r^2 + Mh^2*Wi^2*ddthr*r^2 + 2*Mf^2*Wl^2*ddthl*r^2 - 2*Mf^2*Wl^2*ddthr*r^2 + 4*Mw^2*Wi^2*ddthl*r^2 + 4*Mw^2*Wi^2*ddthr*r^2 - 2*Mf^2*Wl^2*ddthl*r^2*cos(2*thf) + 2*Mf^2*Wl^2*ddthr*r^2*cos(2*thf) - 2*Lb^2*Mb^2*ddthl*r^2*cos(2*thf + 2*thm) + 2*Lb^2*Mb^2*ddthr*r^2*cos(2*thf + 2*thm) + 2*Mb*Mf*Wi^2*ddthl*r^2 + 2*Mb*Mf*Wi^2*ddthr*r^2 + 2*Mb*Mh*Wi^2*ddthl*r^2 + 2*Mb*Mh*Wi^2*ddthr*r^2 + 2*Mf*Mh*Wi^2*ddthl*r^2 + 2*Mf*Mh*Wi^2*ddthr*r^2 + 4*Mb*Mw*Wi^2*ddthl*r^2 + 4*Mb*Mw*Wi^2*ddthr*r^2 + 4*Mf*Mw*Wi^2*ddthl*r^2 + 4*Mf*Mw*Wi^2*ddthr*r^2 + 4*Mh*Mw*Wi^2*ddthl*r^2 + 4*Mh*Mw*Wi^2*ddthr*r^2 - 2*Lb*Mb^2*Wi^2*dthf^2*r*sin(thf + thm) - 2*Lb*Mb^2*Wi^2*dthm^2*r*sin(thf + thm) - 2*Mf^2*Wi^2*Wl*dthf^2*r*sin(thf) + 2*Lb*Mb^2*Wi^2*ddthf*r*cos(thf + thm) + 2*Lb*Mb^2*Wi^2*ddthm*r*cos(thf + thm) + 4*Mf^2*Wl^2*dthf*dthl*r^2*sin(2*thf) - 4*Mf^2*Wl^2*dthf*dthr*r^2*sin(2*thf) + 2*Mf^2*Wi^2*Wl*ddthf*r*cos(thf) + 4*Lb^2*Mb^2*dthf*dthl*r^2*sin(2*thf + 2*thm) - 4*Lb^2*Mb^2*dthf*dthr*r^2*sin(2*thf + 2*thm) + 4*Lb^2*Mb^2*dthl*dthm*r^2*sin(2*thf + 2*thm) - 4*Lb^2*Mb^2*dthm*dthr*r^2*sin(2*thf + 2*thm) - 4*Lb*Mb^2*Wi^2*dthf*dthm*r*sin(thf + thm) - 2*Mb*Mf*Wi^2*Wl*dthf^2*r*sin(thf) - 2*Mf*Mh*Wi^2*Wl*dthf^2*r*sin(thf) - 4*Mf*Mw*Wi^2*Wl*dthf^2*r*sin(thf) + 2*Lb*Mb*Mf*Wi^2*ddthf*r*cos(thf + thm) + 2*Lb*Mb*Mf*Wi^2*ddthm*r*cos(thf + thm) + 2*Lb*Mb*Mh*Wi^2*ddthf*r*cos(thf + thm) + 2*Lb*Mb*Mh*Wi^2*ddthm*r*cos(thf + thm) + 4*Lb*Mb*Mw*Wi^2*ddthf*r*cos(thf + thm) + 4*Lb*Mb*Mw*Wi^2*ddthm*r*cos(thf + thm) + 4*Lb*Mb*Mf*Wl*ddthl*r^2*cos(thm) - 4*Lb*Mb*Mf*Wl*ddthr*r^2*cos(thm) + 2*Mb*Mf*Wi^2*Wl*ddthf*r*cos(thf) + 2*Mf*Mh*Wi^2*Wl*ddthf*r*cos(thf) + 4*Mf*Mw*Wi^2*Wl*ddthf*r*cos(thf) - 4*Lb*Mb*Mf*Wl*ddthl*r^2*cos(2*thf + thm) + 4*Lb*Mb*Mf*Wl*ddthr*r^2*cos(2*thf + thm) - 2*Lb*Mb*Mf*Wi^2*dthf^2*r*sin(thf + thm) - 2*Lb*Mb*Mf*Wi^2*dthm^2*r*sin(thf + thm) - 2*Lb*Mb*Mh*Wi^2*dthf^2*r*sin(thf + thm) - 2*Lb*Mb*Mh*Wi^2*dthm^2*r*sin(thf + thm) - 4*Lb*Mb*Mw*Wi^2*dthf^2*r*sin(thf + thm) - 4*Lb*Mb*Mw*Wi^2*dthm^2*r*sin(thf + thm) - 4*Lb*Mb*Mf*Wl*dthl*dthm*r^2*sin(thm) + 4*Lb*Mb*Mf*Wl*dthm*dthr*r^2*sin(thm) + 8*Lb*Mb*Mf*Wl*dthf*dthl*r^2*sin(2*thf + thm) - 8*Lb*Mb*Mf*Wl*dthf*dthr*r^2*sin(2*thf + thm) + 4*Lb*Mb*Mf*Wl*dthl*dthm*r^2*sin(2*thf + thm) - 4*Lb*Mb*Mf*Wl*dthm*dthr*r^2*sin(2*thf + thm) - 4*Lb*Mb*Mf*Wi^2*dthf*dthm*r*sin(thf + thm) - 4*Lb*Mb*Mh*Wi^2*dthf*dthm*r*sin(thf + thm) - 8*Lb*Mb*Mw*Wi^2*dthf*dthm*r*sin(thf + thm))/(4*Wi^2*(Mb + Mf + Mh + 2*Mw)) == 0;
    (2*Jf*Mb*Wi^2*ddthf + 2*Jf*Mf*Wi^2*ddthf + 2*Jf*Mh*Wi^2*ddthf + 4*Jf*Mw*Wi^2*ddthf + 2*Lb^2*Mb^2*Wi^2*ddthf + 2*Lb^2*Mb^2*Wi^2*ddthm + 2*Mf^2*Wi^2*Wl^2*ddthf - 2*Lb*Mb^2*Wi^2*gravity*sin(thf + thm) - 2*Mf^2*Wi^2*Wl*gravity*sin(thf) - Mf^2*Wl^2*dthl^2*r^2*sin(2*thf) - Mf^2*Wl^2*dthr^2*r^2*sin(2*thf) - Lb^2*Mb^2*dthl^2*r^2*sin(2*thf + 2*thm) - Lb^2*Mb^2*dthr^2*r^2*sin(2*thf + 2*thm) - 2*Lb*Mb*Mf*Wi^2*gravity*sin(thf + thm) - 2*Lb*Mb*Mh*Wi^2*gravity*sin(thf + thm) - 4*Lb*Mb*Mw*Wi^2*gravity*sin(thf + thm) - 2*Mb*Mf*Wi^2*Wl*gravity*sin(thf) - 2*Mf*Mh*Wi^2*Wl*gravity*sin(thf) - 4*Mf*Mw*Wi^2*Wl*gravity*sin(thf) + Lb*Mb^2*Wi^2*ddthl*r*cos(thf + thm) + Lb*Mb^2*Wi^2*ddthr*r*cos(thf + thm) + 2*Mf^2*Wl^2*dthl*dthr*r^2*sin(2*thf) + Mf^2*Wi^2*Wl*ddthl*r*cos(thf) + Mf^2*Wi^2*Wl*ddthr*r*cos(thf) + 2*Lb^2*Mb^2*dthl*dthr*r^2*sin(2*thf + 2*thm) - 2*Lb*Mb*Mf*Wi^2*Wl*dthm^2*sin(thm) - 2*Lb*Mb*Mf*Wl*dthl^2*r^2*sin(2*thf + thm) - 2*Lb*Mb*Mf*Wl*dthr^2*r^2*sin(2*thf + thm) + Lb*Mb*Mf*Wi^2*ddthl*r*cos(thf + thm) + Lb*Mb*Mf*Wi^2*ddthr*r*cos(thf + thm) + Lb*Mb*Mh*Wi^2*ddthl*r*cos(thf + thm) + Lb*Mb*Mh*Wi^2*ddthr*r*cos(thf + thm) + 2*Lb*Mb*Mw*Wi^2*ddthl*r*cos(thf + thm) + 2*Lb*Mb*Mw*Wi^2*ddthr*r*cos(thf + thm) + 4*Lb*Mb*Mf*Wi^2*Wl*ddthf*cos(thm) + 2*Lb*Mb*Mf*Wi^2*Wl*ddthm*cos(thm) + Mb*Mf*Wi^2*Wl*ddthl*r*cos(thf) + Mb*Mf*Wi^2*Wl*ddthr*r*cos(thf) + Mf*Mh*Wi^2*Wl*ddthl*r*cos(thf) + Mf*Mh*Wi^2*Wl*ddthr*r*cos(thf) + 2*Mf*Mw*Wi^2*Wl*ddthl*r*cos(thf) + 2*Mf*Mw*Wi^2*Wl*ddthr*r*cos(thf) - 4*Lb*Mb*Mf*Wi^2*Wl*dthf*dthm*sin(thm) + 4*Lb*Mb*Mf*Wl*dthl*dthr*r^2*sin(2*thf + thm))/(2*Wi^2*(Mb + Mf + Mh + 2*Mw)) == 0
    ];

%It's going to take a lot of time.
sol = solve(eq, ddq);

Answer 1

[Walter Roberson]

Your equations all equal zero, and none of the denominators involve variables that you are solving for. Get rid of the denominators (which is equivalent to multiplying both sides of A/B=0 by B to get A=0)

Now what you have is a series of four equations. Those equations are all linear with respect to the variables you are solving for. You can use https://www.mathworks.com/help/symbolic/sym.equationstomatrix.html and then use \ to get the solution.

[Hlab]

Thank you very much for your quick and excellent answer.

However, when I tried to do it the way you told me to, it took 627 seconds to calculate!

Strangely enough, when I ran it once and then ran the same code again, it took less than a second to calculate.

If you have any other suggestions, please let me know.

clear; format compact
tic
syms Mb Mw Mh Mf real
syms Jh Jw Jb Jh Jf Lb Wi Wl r gravity real
syms thf dthf ddthf real
syms thm dthm ddthm real
syms thr dthr ddthr real
syms thl dthl ddthl real
syms taur taul taub real
syms Jm n Bm real

%the generalized coordinate vector of acceleration
ddq = [ ddthm ddthr ddthl ddthf]';

%equation of motion
eq = [
    (2*Jb*Mb*Wi^2*ddthm - 2*Mf*Wi^2*taub - 2*Mh*Wi^2*taub - 4*Mw*Wi^2*taub - 2*Mb*Wi^2*taub + 2*Jb*Mf*Wi^2*ddthm + 2*Jb*Mh*Wi^2*ddthm + 4*Jb*Mw*Wi^2*ddthm + 2*Lb^2*Mb^2*Wi^2*ddthf + 2*Lb^2*Mb^2*Wi^2*ddthm - 2*Lb*Mb^2*Wi^2*gravity*sin(thf + thm) - Lb^2*Mb^2*dthl^2*r^2*sin(2*thf + 2*thm) - Lb^2*Mb^2*dthr^2*r^2*sin(2*thf + 2*thm) + 2*Bm*Mb*Wi^2*dthm*n^2 + 2*Bm*Mf*Wi^2*dthm*n^2 + 2*Bm*Mh*Wi^2*dthm*n^2 + 4*Bm*Mw*Wi^2*dthm*n^2 - 2*Lb*Mb*Mf*Wi^2*gravity*sin(thf + thm) - 2*Lb*Mb*Mh*Wi^2*gravity*sin(thf + thm) - 4*Lb*Mb*Mw*Wi^2*gravity*sin(thf + thm) + Lb*Mb^2*Wi^2*ddthl*r*cos(thf + thm) + Lb*Mb^2*Wi^2*ddthr*r*cos(thf + thm) + 2*Lb^2*Mb^2*dthl*dthr*r^2*sin(2*thf + 2*thm) + 2*Lb*Mb*Mf*Wi^2*Wl*dthf^2*sin(thm) + Lb*Mb*Mf*Wl*dthl^2*r^2*sin(thm) + Lb*Mb*Mf*Wl*dthr^2*r^2*sin(thm) - Lb*Mb*Mf*Wl*dthl^2*r^2*sin(2*thf + thm) - Lb*Mb*Mf*Wl*dthr^2*r^2*sin(2*thf + thm) + Lb*Mb*Mf*Wi^2*ddthl*r*cos(thf + thm) + Lb*Mb*Mf*Wi^2*ddthr*r*cos(thf + thm) + Lb*Mb*Mh*Wi^2*ddthl*r*cos(thf + thm) + Lb*Mb*Mh*Wi^2*ddthr*r*cos(thf + thm) + 2*Lb*Mb*Mw*Wi^2*ddthl*r*cos(thf + thm) + 2*Lb*Mb*Mw*Wi^2*ddthr*r*cos(thf + thm) + 2*Lb*Mb*Mf*Wi^2*Wl*ddthf*cos(thm) - 2*Lb*Mb*Mf*Wl*dthl*dthr*r^2*sin(thm) + 2*Lb*Mb*Mf*Wl*dthl*dthr*r^2*sin(2*thf + thm)) == 0;
    (4*Jw*Mb*Wi^2*ddthr - 4*Mf*Wi^2*taur - 4*Mh*Wi^2*taur - 8*Mw*Wi^2*taur - 4*Mb*Wi^2*taur + 4*Jw*Mf*Wi^2*ddthr + 4*Jw*Mh*Wi^2*ddthr + 8*Jw*Mw*Wi^2*ddthr - 4*Jh*Mb*ddthl*r^2 + 4*Jh*Mb*ddthr*r^2 - 4*Jh*Mf*ddthl*r^2 + 4*Jh*Mf*ddthr*r^2 - 4*Jh*Mh*ddthl*r^2 + 4*Jh*Mh*ddthr*r^2 - 8*Jh*Mw*ddthl*r^2 + 8*Jh*Mw*ddthr*r^2 - 2*Lb^2*Mb^2*ddthl*r^2 + 2*Lb^2*Mb^2*ddthr*r^2 + Mb^2*Wi^2*ddthl*r^2 + Mb^2*Wi^2*ddthr*r^2 + Mf^2*Wi^2*ddthl*r^2 + Mf^2*Wi^2*ddthr*r^2 + Mh^2*Wi^2*ddthl*r^2 + Mh^2*Wi^2*ddthr*r^2 - 2*Mf^2*Wl^2*ddthl*r^2 + 2*Mf^2*Wl^2*ddthr*r^2 + 4*Mw^2*Wi^2*ddthl*r^2 + 4*Mw^2*Wi^2*ddthr*r^2 + 2*Mf^2*Wl^2*ddthl*r^2*cos(2*thf) - 2*Mf^2*Wl^2*ddthr*r^2*cos(2*thf) + 2*Lb^2*Mb^2*ddthl*r^2*cos(2*thf + 2*thm) - 2*Lb^2*Mb^2*ddthr*r^2*cos(2*thf + 2*thm) + 2*Mb*Mf*Wi^2*ddthl*r^2 + 2*Mb*Mf*Wi^2*ddthr*r^2 + 2*Mb*Mh*Wi^2*ddthl*r^2 + 2*Mb*Mh*Wi^2*ddthr*r^2 + 2*Mf*Mh*Wi^2*ddthl*r^2 + 2*Mf*Mh*Wi^2*ddthr*r^2 + 4*Mb*Mw*Wi^2*ddthl*r^2 + 4*Mb*Mw*Wi^2*ddthr*r^2 + 4*Mf*Mw*Wi^2*ddthl*r^2 + 4*Mf*Mw*Wi^2*ddthr*r^2 + 4*Mh*Mw*Wi^2*ddthl*r^2 + 4*Mh*Mw*Wi^2*ddthr*r^2 - 2*Lb*Mb^2*Wi^2*dthf^2*r*sin(thf + thm) - 2*Lb*Mb^2*Wi^2*dthm^2*r*sin(thf + thm) - 2*Mf^2*Wi^2*Wl*dthf^2*r*sin(thf) + 2*Lb*Mb^2*Wi^2*ddthf*r*cos(thf + thm) + 2*Lb*Mb^2*Wi^2*ddthm*r*cos(thf + thm) - 4*Mf^2*Wl^2*dthf*dthl*r^2*sin(2*thf) + 4*Mf^2*Wl^2*dthf*dthr*r^2*sin(2*thf) + 2*Mf^2*Wi^2*Wl*ddthf*r*cos(thf) - 4*Lb^2*Mb^2*dthf*dthl*r^2*sin(2*thf + 2*thm) + 4*Lb^2*Mb^2*dthf*dthr*r^2*sin(2*thf + 2*thm) - 4*Lb^2*Mb^2*dthl*dthm*r^2*sin(2*thf + 2*thm) + 4*Lb^2*Mb^2*dthm*dthr*r^2*sin(2*thf + 2*thm) - 4*Lb*Mb^2*Wi^2*dthf*dthm*r*sin(thf + thm) - 2*Mb*Mf*Wi^2*Wl*dthf^2*r*sin(thf) - 2*Mf*Mh*Wi^2*Wl*dthf^2*r*sin(thf) - 4*Mf*Mw*Wi^2*Wl*dthf^2*r*sin(thf) + 2*Lb*Mb*Mf*Wi^2*ddthf*r*cos(thf + thm) + 2*Lb*Mb*Mf*Wi^2*ddthm*r*cos(thf + thm) + 2*Lb*Mb*Mh*Wi^2*ddthf*r*cos(thf + thm) + 2*Lb*Mb*Mh*Wi^2*ddthm*r*cos(thf + thm) + 4*Lb*Mb*Mw*Wi^2*ddthf*r*cos(thf + thm) + 4*Lb*Mb*Mw*Wi^2*ddthm*r*cos(thf + thm) - 4*Lb*Mb*Mf*Wl*ddthl*r^2*cos(thm) + 4*Lb*Mb*Mf*Wl*ddthr*r^2*cos(thm) + 2*Mb*Mf*Wi^2*Wl*ddthf*r*cos(thf) + 2*Mf*Mh*Wi^2*Wl*ddthf*r*cos(thf) + 4*Mf*Mw*Wi^2*Wl*ddthf*r*cos(thf) + 4*Lb*Mb*Mf*Wl*ddthl*r^2*cos(2*thf + thm) - 4*Lb*Mb*Mf*Wl*ddthr*r^2*cos(2*thf + thm) - 2*Lb*Mb*Mf*Wi^2*dthf^2*r*sin(thf + thm) - 2*Lb*Mb*Mf*Wi^2*dthm^2*r*sin(thf + thm) - 2*Lb*Mb*Mh*Wi^2*dthf^2*r*sin(thf + thm) - 2*Lb*Mb*Mh*Wi^2*dthm^2*r*sin(thf + thm) - 4*Lb*Mb*Mw*Wi^2*dthf^2*r*sin(thf + thm) - 4*Lb*Mb*Mw*Wi^2*dthm^2*r*sin(thf + thm) + 4*Lb*Mb*Mf*Wl*dthl*dthm*r^2*sin(thm) - 4*Lb*Mb*Mf*Wl*dthm*dthr*r^2*sin(thm) - 8*Lb*Mb*Mf*Wl*dthf*dthl*r^2*sin(2*thf + thm) + 8*Lb*Mb*Mf*Wl*dthf*dthr*r^2*sin(2*thf + thm) - 4*Lb*Mb*Mf*Wl*dthl*dthm*r^2*sin(2*thf + thm) + 4*Lb*Mb*Mf*Wl*dthm*dthr*r^2*sin(2*thf + thm) - 4*Lb*Mb*Mf*Wi^2*dthf*dthm*r*sin(thf + thm) - 4*Lb*Mb*Mh*Wi^2*dthf*dthm*r*sin(thf + thm) - 8*Lb*Mb*Mw*Wi^2*dthf*dthm*r*sin(thf + thm)) == 0;
    (4*Jw*Mb*Wi^2*ddthl - 4*Mf*Wi^2*taul - 4*Mh*Wi^2*taul - 8*Mw*Wi^2*taul - 4*Mb*Wi^2*taul + 4*Jw*Mf*Wi^2*ddthl + 4*Jw*Mh*Wi^2*ddthl + 8*Jw*Mw*Wi^2*ddthl + 4*Jh*Mb*ddthl*r^2 - 4*Jh*Mb*ddthr*r^2 + 4*Jh*Mf*ddthl*r^2 - 4*Jh*Mf*ddthr*r^2 + 4*Jh*Mh*ddthl*r^2 - 4*Jh*Mh*ddthr*r^2 + 8*Jh*Mw*ddthl*r^2 - 8*Jh*Mw*ddthr*r^2 + 2*Lb^2*Mb^2*ddthl*r^2 - 2*Lb^2*Mb^2*ddthr*r^2 + Mb^2*Wi^2*ddthl*r^2 + Mb^2*Wi^2*ddthr*r^2 + Mf^2*Wi^2*ddthl*r^2 + Mf^2*Wi^2*ddthr*r^2 + Mh^2*Wi^2*ddthl*r^2 + Mh^2*Wi^2*ddthr*r^2 + 2*Mf^2*Wl^2*ddthl*r^2 - 2*Mf^2*Wl^2*ddthr*r^2 + 4*Mw^2*Wi^2*ddthl*r^2 + 4*Mw^2*Wi^2*ddthr*r^2 - 2*Mf^2*Wl^2*ddthl*r^2*cos(2*thf) + 2*Mf^2*Wl^2*ddthr*r^2*cos(2*thf) - 2*Lb^2*Mb^2*ddthl*r^2*cos(2*thf + 2*thm) + 2*Lb^2*Mb^2*ddthr*r^2*cos(2*thf + 2*thm) + 2*Mb*Mf*Wi^2*ddthl*r^2 + 2*Mb*Mf*Wi^2*ddthr*r^2 + 2*Mb*Mh*Wi^2*ddthl*r^2 + 2*Mb*Mh*Wi^2*ddthr*r^2 + 2*Mf*Mh*Wi^2*ddthl*r^2 + 2*Mf*Mh*Wi^2*ddthr*r^2 + 4*Mb*Mw*Wi^2*ddthl*r^2 + 4*Mb*Mw*Wi^2*ddthr*r^2 + 4*Mf*Mw*Wi^2*ddthl*r^2 + 4*Mf*Mw*Wi^2*ddthr*r^2 + 4*Mh*Mw*Wi^2*ddthl*r^2 + 4*Mh*Mw*Wi^2*ddthr*r^2 - 2*Lb*Mb^2*Wi^2*dthf^2*r*sin(thf + thm) - 2*Lb*Mb^2*Wi^2*dthm^2*r*sin(thf + thm) - 2*Mf^2*Wi^2*Wl*dthf^2*r*sin(thf) + 2*Lb*Mb^2*Wi^2*ddthf*r*cos(thf + thm) + 2*Lb*Mb^2*Wi^2*ddthm*r*cos(thf + thm) + 4*Mf^2*Wl^2*dthf*dthl*r^2*sin(2*thf) - 4*Mf^2*Wl^2*dthf*dthr*r^2*sin(2*thf) + 2*Mf^2*Wi^2*Wl*ddthf*r*cos(thf) + 4*Lb^2*Mb^2*dthf*dthl*r^2*sin(2*thf + 2*thm) - 4*Lb^2*Mb^2*dthf*dthr*r^2*sin(2*thf + 2*thm) + 4*Lb^2*Mb^2*dthl*dthm*r^2*sin(2*thf + 2*thm) - 4*Lb^2*Mb^2*dthm*dthr*r^2*sin(2*thf + 2*thm) - 4*Lb*Mb^2*Wi^2*dthf*dthm*r*sin(thf + thm) - 2*Mb*Mf*Wi^2*Wl*dthf^2*r*sin(thf) - 2*Mf*Mh*Wi^2*Wl*dthf^2*r*sin(thf) - 4*Mf*Mw*Wi^2*Wl*dthf^2*r*sin(thf) + 2*Lb*Mb*Mf*Wi^2*ddthf*r*cos(thf + thm) + 2*Lb*Mb*Mf*Wi^2*ddthm*r*cos(thf + thm) + 2*Lb*Mb*Mh*Wi^2*ddthf*r*cos(thf + thm) + 2*Lb*Mb*Mh*Wi^2*ddthm*r*cos(thf + thm) + 4*Lb*Mb*Mw*Wi^2*ddthf*r*cos(thf + thm) + 4*Lb*Mb*Mw*Wi^2*ddthm*r*cos(thf + thm) + 4*Lb*Mb*Mf*Wl*ddthl*r^2*cos(thm) - 4*Lb*Mb*Mf*Wl*ddthr*r^2*cos(thm) + 2*Mb*Mf*Wi^2*Wl*ddthf*r*cos(thf) + 2*Mf*Mh*Wi^2*Wl*ddthf*r*cos(thf) + 4*Mf*Mw*Wi^2*Wl*ddthf*r*cos(thf) - 4*Lb*Mb*Mf*Wl*ddthl*r^2*cos(2*thf + thm) + 4*Lb*Mb*Mf*Wl*ddthr*r^2*cos(2*thf + thm) - 2*Lb*Mb*Mf*Wi^2*dthf^2*r*sin(thf + thm) - 2*Lb*Mb*Mf*Wi^2*dthm^2*r*sin(thf + thm) - 2*Lb*Mb*Mh*Wi^2*dthf^2*r*sin(thf + thm) - 2*Lb*Mb*Mh*Wi^2*dthm^2*r*sin(thf + thm) - 4*Lb*Mb*Mw*Wi^2*dthf^2*r*sin(thf + thm) - 4*Lb*Mb*Mw*Wi^2*dthm^2*r*sin(thf + thm) - 4*Lb*Mb*Mf*Wl*dthl*dthm*r^2*sin(thm) + 4*Lb*Mb*Mf*Wl*dthm*dthr*r^2*sin(thm) + 8*Lb*Mb*Mf*Wl*dthf*dthl*r^2*sin(2*thf + thm) - 8*Lb*Mb*Mf*Wl*dthf*dthr*r^2*sin(2*thf + thm) + 4*Lb*Mb*Mf*Wl*dthl*dthm*r^2*sin(2*thf + thm) - 4*Lb*Mb*Mf*Wl*dthm*dthr*r^2*sin(2*thf + thm) - 4*Lb*Mb*Mf*Wi^2*dthf*dthm*r*sin(thf + thm) - 4*Lb*Mb*Mh*Wi^2*dthf*dthm*r*sin(thf + thm) - 8*Lb*Mb*Mw*Wi^2*dthf*dthm*r*sin(thf + thm)) == 0;
    (2*Jf*Mb*Wi^2*ddthf + 2*Jf*Mf*Wi^2*ddthf + 2*Jf*Mh*Wi^2*ddthf + 4*Jf*Mw*Wi^2*ddthf + 2*Lb^2*Mb^2*Wi^2*ddthf + 2*Lb^2*Mb^2*Wi^2*ddthm + 2*Mf^2*Wi^2*Wl^2*ddthf - 2*Lb*Mb^2*Wi^2*gravity*sin(thf + thm) - 2*Mf^2*Wi^2*Wl*gravity*sin(thf) - Mf^2*Wl^2*dthl^2*r^2*sin(2*thf) - Mf^2*Wl^2*dthr^2*r^2*sin(2*thf) - Lb^2*Mb^2*dthl^2*r^2*sin(2*thf + 2*thm) - Lb^2*Mb^2*dthr^2*r^2*sin(2*thf + 2*thm) - 2*Lb*Mb*Mf*Wi^2*gravity*sin(thf + thm) - 2*Lb*Mb*Mh*Wi^2*gravity*sin(thf + thm) - 4*Lb*Mb*Mw*Wi^2*gravity*sin(thf + thm) - 2*Mb*Mf*Wi^2*Wl*gravity*sin(thf) - 2*Mf*Mh*Wi^2*Wl*gravity*sin(thf) - 4*Mf*Mw*Wi^2*Wl*gravity*sin(thf) + Lb*Mb^2*Wi^2*ddthl*r*cos(thf + thm) + Lb*Mb^2*Wi^2*ddthr*r*cos(thf + thm) + 2*Mf^2*Wl^2*dthl*dthr*r^2*sin(2*thf) + Mf^2*Wi^2*Wl*ddthl*r*cos(thf) + Mf^2*Wi^2*Wl*ddthr*r*cos(thf) + 2*Lb^2*Mb^2*dthl*dthr*r^2*sin(2*thf + 2*thm) - 2*Lb*Mb*Mf*Wi^2*Wl*dthm^2*sin(thm) - 2*Lb*Mb*Mf*Wl*dthl^2*r^2*sin(2*thf + thm) - 2*Lb*Mb*Mf*Wl*dthr^2*r^2*sin(2*thf + thm) + Lb*Mb*Mf*Wi^2*ddthl*r*cos(thf + thm) + Lb*Mb*Mf*Wi^2*ddthr*r*cos(thf + thm) + Lb*Mb*Mh*Wi^2*ddthl*r*cos(thf + thm) + Lb*Mb*Mh*Wi^2*ddthr*r*cos(thf + thm) + 2*Lb*Mb*Mw*Wi^2*ddthl*r*cos(thf + thm) + 2*Lb*Mb*Mw*Wi^2*ddthr*r*cos(thf + thm) + 4*Lb*Mb*Mf*Wi^2*Wl*ddthf*cos(thm) + 2*Lb*Mb*Mf*Wi^2*Wl*ddthm*cos(thm) + Mb*Mf*Wi^2*Wl*ddthl*r*cos(thf) + Mb*Mf*Wi^2*Wl*ddthr*r*cos(thf) + Mf*Mh*Wi^2*Wl*ddthl*r*cos(thf) + Mf*Mh*Wi^2*Wl*ddthr*r*cos(thf) + 2*Mf*Mw*Wi^2*Wl*ddthl*r*cos(thf) + 2*Mf*Mw*Wi^2*Wl*ddthr*r*cos(thf) - 4*Lb*Mb*Mf*Wi^2*Wl*dthf*dthm*sin(thm) + 4*Lb*Mb*Mf*Wl*dthl*dthr*r^2*sin(2*thf + thm)) == 0
    ];

[A,b] = equationsToMatrix(eq,ddq);
sol = A\b;
toc
% 
% %It's going to take a lot of time.
% sol = solve(eq, ddq);
% simplify(sol.ddthm)

[Walter Roberson]

There is a trick.

Create a generic 4 x 4 symbolic matrix

syms As [4 4]
syms bs [4 1]
sols = As\bs;

Now that you have the generic form

sol = subs(sols, {As, bs}, {A, b});

[Hlab]

Thank you for sharing this trick with me.

My calculation time has been reduced from 600 seconds to about 2 seconds.

I've been struggling with this problem for a long time and I'm thrilled to have solved it.

I would also like to know why this happens if possible.

[Walter Roberson]

The general case for solving linear equations is an extended formula,

syms As [4 4]
syms bs [4 1]
sols = As\bs

sols = 

A lot of input variables getting multiplied by other input variables. It is pretty much calculating the determinant and calculating the inverse based on the determinant, and multiplying that by the b vector.

When the terms are messy:

syms Mb Mw Mh Mf real
syms Jh Jw Jb Jh Jf Lb Wi Wl r gravity real
syms thf dthf ddthf real
syms thm dthm ddthm real
syms thr dthr ddthr real
syms thl dthl ddthl real
syms taur taul taub real
syms Jm n Bm real
%the generalized coordinate vector of acceleration
ddq = [ ddthm ddthr ddthl ddthf]';
%equation of motion
eq = [
    (2*Jb*Mb*Wi^2*ddthm - 2*Mf*Wi^2*taub - 2*Mh*Wi^2*taub - 4*Mw*Wi^2*taub - 2*Mb*Wi^2*taub + 2*Jb*Mf*Wi^2*ddthm + 2*Jb*Mh*Wi^2*ddthm + 4*Jb*Mw*Wi^2*ddthm + 2*Lb^2*Mb^2*Wi^2*ddthf + 2*Lb^2*Mb^2*Wi^2*ddthm - 2*Lb*Mb^2*Wi^2*gravity*sin(thf + thm) - Lb^2*Mb^2*dthl^2*r^2*sin(2*thf + 2*thm) - Lb^2*Mb^2*dthr^2*r^2*sin(2*thf + 2*thm) + 2*Bm*Mb*Wi^2*dthm*n^2 + 2*Bm*Mf*Wi^2*dthm*n^2 + 2*Bm*Mh*Wi^2*dthm*n^2 + 4*Bm*Mw*Wi^2*dthm*n^2 - 2*Lb*Mb*Mf*Wi^2*gravity*sin(thf + thm) - 2*Lb*Mb*Mh*Wi^2*gravity*sin(thf + thm) - 4*Lb*Mb*Mw*Wi^2*gravity*sin(thf + thm) + Lb*Mb^2*Wi^2*ddthl*r*cos(thf + thm) + Lb*Mb^2*Wi^2*ddthr*r*cos(thf + thm) + 2*Lb^2*Mb^2*dthl*dthr*r^2*sin(2*thf + 2*thm) + 2*Lb*Mb*Mf*Wi^2*Wl*dthf^2*sin(thm) + Lb*Mb*Mf*Wl*dthl^2*r^2*sin(thm) + Lb*Mb*Mf*Wl*dthr^2*r^2*sin(thm) - Lb*Mb*Mf*Wl*dthl^2*r^2*sin(2*thf + thm) - Lb*Mb*Mf*Wl*dthr^2*r^2*sin(2*thf + thm) + Lb*Mb*Mf*Wi^2*ddthl*r*cos(thf + thm) + Lb*Mb*Mf*Wi^2*ddthr*r*cos(thf + thm) + Lb*Mb*Mh*Wi^2*ddthl*r*cos(thf + thm) + Lb*Mb*Mh*Wi^2*ddthr*r*cos(thf + thm) + 2*Lb*Mb*Mw*Wi^2*ddthl*r*cos(thf + thm) + 2*Lb*Mb*Mw*Wi^2*ddthr*r*cos(thf + thm) + 2*Lb*Mb*Mf*Wi^2*Wl*ddthf*cos(thm) - 2*Lb*Mb*Mf*Wl*dthl*dthr*r^2*sin(thm) + 2*Lb*Mb*Mf*Wl*dthl*dthr*r^2*sin(2*thf + thm)) == 0;
    (4*Jw*Mb*Wi^2*ddthr - 4*Mf*Wi^2*taur - 4*Mh*Wi^2*taur - 8*Mw*Wi^2*taur - 4*Mb*Wi^2*taur + 4*Jw*Mf*Wi^2*ddthr + 4*Jw*Mh*Wi^2*ddthr + 8*Jw*Mw*Wi^2*ddthr - 4*Jh*Mb*ddthl*r^2 + 4*Jh*Mb*ddthr*r^2 - 4*Jh*Mf*ddthl*r^2 + 4*Jh*Mf*ddthr*r^2 - 4*Jh*Mh*ddthl*r^2 + 4*Jh*Mh*ddthr*r^2 - 8*Jh*Mw*ddthl*r^2 + 8*Jh*Mw*ddthr*r^2 - 2*Lb^2*Mb^2*ddthl*r^2 + 2*Lb^2*Mb^2*ddthr*r^2 + Mb^2*Wi^2*ddthl*r^2 + Mb^2*Wi^2*ddthr*r^2 + Mf^2*Wi^2*ddthl*r^2 + Mf^2*Wi^2*ddthr*r^2 + Mh^2*Wi^2*ddthl*r^2 + Mh^2*Wi^2*ddthr*r^2 - 2*Mf^2*Wl^2*ddthl*r^2 + 2*Mf^2*Wl^2*ddthr*r^2 + 4*Mw^2*Wi^2*ddthl*r^2 + 4*Mw^2*Wi^2*ddthr*r^2 + 2*Mf^2*Wl^2*ddthl*r^2*cos(2*thf) - 2*Mf^2*Wl^2*ddthr*r^2*cos(2*thf) + 2*Lb^2*Mb^2*ddthl*r^2*cos(2*thf + 2*thm) - 2*Lb^2*Mb^2*ddthr*r^2*cos(2*thf + 2*thm) + 2*Mb*Mf*Wi^2*ddthl*r^2 + 2*Mb*Mf*Wi^2*ddthr*r^2 + 2*Mb*Mh*Wi^2*ddthl*r^2 + 2*Mb*Mh*Wi^2*ddthr*r^2 + 2*Mf*Mh*Wi^2*ddthl*r^2 + 2*Mf*Mh*Wi^2*ddthr*r^2 + 4*Mb*Mw*Wi^2*ddthl*r^2 + 4*Mb*Mw*Wi^2*ddthr*r^2 + 4*Mf*Mw*Wi^2*ddthl*r^2 + 4*Mf*Mw*Wi^2*ddthr*r^2 + 4*Mh*Mw*Wi^2*ddthl*r^2 + 4*Mh*Mw*Wi^2*ddthr*r^2 - 2*Lb*Mb^2*Wi^2*dthf^2*r*sin(thf + thm) - 2*Lb*Mb^2*Wi^2*dthm^2*r*sin(thf + thm) - 2*Mf^2*Wi^2*Wl*dthf^2*r*sin(thf) + 2*Lb*Mb^2*Wi^2*ddthf*r*cos(thf + thm) + 2*Lb*Mb^2*Wi^2*ddthm*r*cos(thf + thm) - 4*Mf^2*Wl^2*dthf*dthl*r^2*sin(2*thf) + 4*Mf^2*Wl^2*dthf*dthr*r^2*sin(2*thf) + 2*Mf^2*Wi^2*Wl*ddthf*r*cos(thf) - 4*Lb^2*Mb^2*dthf*dthl*r^2*sin(2*thf + 2*thm) + 4*Lb^2*Mb^2*dthf*dthr*r^2*sin(2*thf + 2*thm) - 4*Lb^2*Mb^2*dthl*dthm*r^2*sin(2*thf + 2*thm) + 4*Lb^2*Mb^2*dthm*dthr*r^2*sin(2*thf + 2*thm) - 4*Lb*Mb^2*Wi^2*dthf*dthm*r*sin(thf + thm) - 2*Mb*Mf*Wi^2*Wl*dthf^2*r*sin(thf) - 2*Mf*Mh*Wi^2*Wl*dthf^2*r*sin(thf) - 4*Mf*Mw*Wi^2*Wl*dthf^2*r*sin(thf) + 2*Lb*Mb*Mf*Wi^2*ddthf*r*cos(thf + thm) + 2*Lb*Mb*Mf*Wi^2*ddthm*r*cos(thf + thm) + 2*Lb*Mb*Mh*Wi^2*ddthf*r*cos(thf + thm) + 2*Lb*Mb*Mh*Wi^2*ddthm*r*cos(thf + thm) + 4*Lb*Mb*Mw*Wi^2*ddthf*r*cos(thf + thm) + 4*Lb*Mb*Mw*Wi^2*ddthm*r*cos(thf + thm) - 4*Lb*Mb*Mf*Wl*ddthl*r^2*cos(thm) + 4*Lb*Mb*Mf*Wl*ddthr*r^2*cos(thm) + 2*Mb*Mf*Wi^2*Wl*ddthf*r*cos(thf) + 2*Mf*Mh*Wi^2*Wl*ddthf*r*cos(thf) + 4*Mf*Mw*Wi^2*Wl*ddthf*r*cos(thf) + 4*Lb*Mb*Mf*Wl*ddthl*r^2*cos(2*thf + thm) - 4*Lb*Mb*Mf*Wl*ddthr*r^2*cos(2*thf + thm) - 2*Lb*Mb*Mf*Wi^2*dthf^2*r*sin(thf + thm) - 2*Lb*Mb*Mf*Wi^2*dthm^2*r*sin(thf + thm) - 2*Lb*Mb*Mh*Wi^2*dthf^2*r*sin(thf + thm) - 2*Lb*Mb*Mh*Wi^2*dthm^2*r*sin(thf + thm) - 4*Lb*Mb*Mw*Wi^2*dthf^2*r*sin(thf + thm) - 4*Lb*Mb*Mw*Wi^2*dthm^2*r*sin(thf + thm) + 4*Lb*Mb*Mf*Wl*dthl*dthm*r^2*sin(thm) - 4*Lb*Mb*Mf*Wl*dthm*dthr*r^2*sin(thm) - 8*Lb*Mb*Mf*Wl*dthf*dthl*r^2*sin(2*thf + thm) + 8*Lb*Mb*Mf*Wl*dthf*dthr*r^2*sin(2*thf + thm) - 4*Lb*Mb*Mf*Wl*dthl*dthm*r^2*sin(2*thf + thm) + 4*Lb*Mb*Mf*Wl*dthm*dthr*r^2*sin(2*thf + thm) - 4*Lb*Mb*Mf*Wi^2*dthf*dthm*r*sin(thf + thm) - 4*Lb*Mb*Mh*Wi^2*dthf*dthm*r*sin(thf + thm) - 8*Lb*Mb*Mw*Wi^2*dthf*dthm*r*sin(thf + thm)) == 0;
    (4*Jw*Mb*Wi^2*ddthl - 4*Mf*Wi^2*taul - 4*Mh*Wi^2*taul - 8*Mw*Wi^2*taul - 4*Mb*Wi^2*taul + 4*Jw*Mf*Wi^2*ddthl + 4*Jw*Mh*Wi^2*ddthl + 8*Jw*Mw*Wi^2*ddthl + 4*Jh*Mb*ddthl*r^2 - 4*Jh*Mb*ddthr*r^2 + 4*Jh*Mf*ddthl*r^2 - 4*Jh*Mf*ddthr*r^2 + 4*Jh*Mh*ddthl*r^2 - 4*Jh*Mh*ddthr*r^2 + 8*Jh*Mw*ddthl*r^2 - 8*Jh*Mw*ddthr*r^2 + 2*Lb^2*Mb^2*ddthl*r^2 - 2*Lb^2*Mb^2*ddthr*r^2 + Mb^2*Wi^2*ddthl*r^2 + Mb^2*Wi^2*ddthr*r^2 + Mf^2*Wi^2*ddthl*r^2 + Mf^2*Wi^2*ddthr*r^2 + Mh^2*Wi^2*ddthl*r^2 + Mh^2*Wi^2*ddthr*r^2 + 2*Mf^2*Wl^2*ddthl*r^2 - 2*Mf^2*Wl^2*ddthr*r^2 + 4*Mw^2*Wi^2*ddthl*r^2 + 4*Mw^2*Wi^2*ddthr*r^2 - 2*Mf^2*Wl^2*ddthl*r^2*cos(2*thf) + 2*Mf^2*Wl^2*ddthr*r^2*cos(2*thf) - 2*Lb^2*Mb^2*ddthl*r^2*cos(2*thf + 2*thm) + 2*Lb^2*Mb^2*ddthr*r^2*cos(2*thf + 2*thm) + 2*Mb*Mf*Wi^2*ddthl*r^2 + 2*Mb*Mf*Wi^2*ddthr*r^2 + 2*Mb*Mh*Wi^2*ddthl*r^2 + 2*Mb*Mh*Wi^2*ddthr*r^2 + 2*Mf*Mh*Wi^2*ddthl*r^2 + 2*Mf*Mh*Wi^2*ddthr*r^2 + 4*Mb*Mw*Wi^2*ddthl*r^2 + 4*Mb*Mw*Wi^2*ddthr*r^2 + 4*Mf*Mw*Wi^2*ddthl*r^2 + 4*Mf*Mw*Wi^2*ddthr*r^2 + 4*Mh*Mw*Wi^2*ddthl*r^2 + 4*Mh*Mw*Wi^2*ddthr*r^2 - 2*Lb*Mb^2*Wi^2*dthf^2*r*sin(thf + thm) - 2*Lb*Mb^2*Wi^2*dthm^2*r*sin(thf + thm) - 2*Mf^2*Wi^2*Wl*dthf^2*r*sin(thf) + 2*Lb*Mb^2*Wi^2*ddthf*r*cos(thf + thm) + 2*Lb*Mb^2*Wi^2*ddthm*r*cos(thf + thm) + 4*Mf^2*Wl^2*dthf*dthl*r^2*sin(2*thf) - 4*Mf^2*Wl^2*dthf*dthr*r^2*sin(2*thf) + 2*Mf^2*Wi^2*Wl*ddthf*r*cos(thf) + 4*Lb^2*Mb^2*dthf*dthl*r^2*sin(2*thf + 2*thm) - 4*Lb^2*Mb^2*dthf*dthr*r^2*sin(2*thf + 2*thm) + 4*Lb^2*Mb^2*dthl*dthm*r^2*sin(2*thf + 2*thm) - 4*Lb^2*Mb^2*dthm*dthr*r^2*sin(2*thf + 2*thm) - 4*Lb*Mb^2*Wi^2*dthf*dthm*r*sin(thf + thm) - 2*Mb*Mf*Wi^2*Wl*dthf^2*r*sin(thf) - 2*Mf*Mh*Wi^2*Wl*dthf^2*r*sin(thf) - 4*Mf*Mw*Wi^2*Wl*dthf^2*r*sin(thf) + 2*Lb*Mb*Mf*Wi^2*ddthf*r*cos(thf + thm) + 2*Lb*Mb*Mf*Wi^2*ddthm*r*cos(thf + thm) + 2*Lb*Mb*Mh*Wi^2*ddthf*r*cos(thf + thm) + 2*Lb*Mb*Mh*Wi^2*ddthm*r*cos(thf + thm) + 4*Lb*Mb*Mw*Wi^2*ddthf*r*cos(thf + thm) + 4*Lb*Mb*Mw*Wi^2*ddthm*r*cos(thf + thm) + 4*Lb*Mb*Mf*Wl*ddthl*r^2*cos(thm) - 4*Lb*Mb*Mf*Wl*ddthr*r^2*cos(thm) + 2*Mb*Mf*Wi^2*Wl*ddthf*r*cos(thf) + 2*Mf*Mh*Wi^2*Wl*ddthf*r*cos(thf) + 4*Mf*Mw*Wi^2*Wl*ddthf*r*cos(thf) - 4*Lb*Mb*Mf*Wl*ddthl*r^2*cos(2*thf + thm) + 4*Lb*Mb*Mf*Wl*ddthr*r^2*cos(2*thf + thm) - 2*Lb*Mb*Mf*Wi^2*dthf^2*r*sin(thf + thm) - 2*Lb*Mb*Mf*Wi^2*dthm^2*r*sin(thf + thm) - 2*Lb*Mb*Mh*Wi^2*dthf^2*r*sin(thf + thm) - 2*Lb*Mb*Mh*Wi^2*dthm^2*r*sin(thf + thm) - 4*Lb*Mb*Mw*Wi^2*dthf^2*r*sin(thf + thm) - 4*Lb*Mb*Mw*Wi^2*dthm^2*r*sin(thf + thm) - 4*Lb*Mb*Mf*Wl*dthl*dthm*r^2*sin(thm) + 4*Lb*Mb*Mf*Wl*dthm*dthr*r^2*sin(thm) + 8*Lb*Mb*Mf*Wl*dthf*dthl*r^2*sin(2*thf + thm) - 8*Lb*Mb*Mf*Wl*dthf*dthr*r^2*sin(2*thf + thm) + 4*Lb*Mb*Mf*Wl*dthl*dthm*r^2*sin(2*thf + thm) - 4*Lb*Mb*Mf*Wl*dthm*dthr*r^2*sin(2*thf + thm) - 4*Lb*Mb*Mf*Wi^2*dthf*dthm*r*sin(thf + thm) - 4*Lb*Mb*Mh*Wi^2*dthf*dthm*r*sin(thf + thm) - 8*Lb*Mb*Mw*Wi^2*dthf*dthm*r*sin(thf + thm)) == 0;
    (2*Jf*Mb*Wi^2*ddthf + 2*Jf*Mf*Wi^2*ddthf + 2*Jf*Mh*Wi^2*ddthf + 4*Jf*Mw*Wi^2*ddthf + 2*Lb^2*Mb^2*Wi^2*ddthf + 2*Lb^2*Mb^2*Wi^2*ddthm + 2*Mf^2*Wi^2*Wl^2*ddthf - 2*Lb*Mb^2*Wi^2*gravity*sin(thf + thm) - 2*Mf^2*Wi^2*Wl*gravity*sin(thf) - Mf^2*Wl^2*dthl^2*r^2*sin(2*thf) - Mf^2*Wl^2*dthr^2*r^2*sin(2*thf) - Lb^2*Mb^2*dthl^2*r^2*sin(2*thf + 2*thm) - Lb^2*Mb^2*dthr^2*r^2*sin(2*thf + 2*thm) - 2*Lb*Mb*Mf*Wi^2*gravity*sin(thf + thm) - 2*Lb*Mb*Mh*Wi^2*gravity*sin(thf + thm) - 4*Lb*Mb*Mw*Wi^2*gravity*sin(thf + thm) - 2*Mb*Mf*Wi^2*Wl*gravity*sin(thf) - 2*Mf*Mh*Wi^2*Wl*gravity*sin(thf) - 4*Mf*Mw*Wi^2*Wl*gravity*sin(thf) + Lb*Mb^2*Wi^2*ddthl*r*cos(thf + thm) + Lb*Mb^2*Wi^2*ddthr*r*cos(thf + thm) + 2*Mf^2*Wl^2*dthl*dthr*r^2*sin(2*thf) + Mf^2*Wi^2*Wl*ddthl*r*cos(thf) + Mf^2*Wi^2*Wl*ddthr*r*cos(thf) + 2*Lb^2*Mb^2*dthl*dthr*r^2*sin(2*thf + 2*thm) - 2*Lb*Mb*Mf*Wi^2*Wl*dthm^2*sin(thm) - 2*Lb*Mb*Mf*Wl*dthl^2*r^2*sin(2*thf + thm) - 2*Lb*Mb*Mf*Wl*dthr^2*r^2*sin(2*thf + thm) + Lb*Mb*Mf*Wi^2*ddthl*r*cos(thf + thm) + Lb*Mb*Mf*Wi^2*ddthr*r*cos(thf + thm) + Lb*Mb*Mh*Wi^2*ddthl*r*cos(thf + thm) + Lb*Mb*Mh*Wi^2*ddthr*r*cos(thf + thm) + 2*Lb*Mb*Mw*Wi^2*ddthl*r*cos(thf + thm) + 2*Lb*Mb*Mw*Wi^2*ddthr*r*cos(thf + thm) + 4*Lb*Mb*Mf*Wi^2*Wl*ddthf*cos(thm) + 2*Lb*Mb*Mf*Wi^2*Wl*ddthm*cos(thm) + Mb*Mf*Wi^2*Wl*ddthl*r*cos(thf) + Mb*Mf*Wi^2*Wl*ddthr*r*cos(thf) + Mf*Mh*Wi^2*Wl*ddthl*r*cos(thf) + Mf*Mh*Wi^2*Wl*ddthr*r*cos(thf) + 2*Mf*Mw*Wi^2*Wl*ddthl*r*cos(thf) + 2*Mf*Mw*Wi^2*Wl*ddthr*r*cos(thf) - 4*Lb*Mb*Mf*Wi^2*Wl*dthf*dthm*sin(thm) + 4*Lb*Mb*Mf*Wl*dthl*dthr*r^2*sin(2*thf + thm)) == 0
    ];
[A,b] = equationsToMatrix(eq,ddq)

A = 

b = 

Then the individual entries get messy, such as

A(1,1)

ans = 

and multiple such terms have to get multiplied through and added to get the results. But as it goes through, the \ operator ends up pretty much having to expand() all of those multiplications, which results in large internal terms, and it has to fish through the expanded versions to figure out what can be added together... which involves comparing every term to every other term in the addition to determine whether they happen to have the same set of variables to the same powers and so should be added together. It is a lot of work to carry around from step to step.

When you compute the generic form with As and bs and subs() in the actual terms, then within any one multiplication it has to collect powers (so A*B*A --> A^2*B) but it does not have to expand out all of the multiplications. (P+Q)*(R+S)*(T+U) gets left in that form instead of being multiplied out like the internal processing of the \ operator had to do in order to be sure it was getting its computation right.