As shown in the following matlab program which is used to find numerical solutions of a dynamic equations (differential equations) which are developed by Lagrange method, the ode45, as well as ode23s, ode15s solvers in matlab are used solve the equations. However, the program has kept busy for even 8 hours. How could I deal with this problem?
the main M file is as follows:
clear; clc; close;
% constant
rouL = 2550;
L1 = 0.2;
L2 = 0.2;
bL = 0.015;
hL = 0.005;
mL = rouL*bL*hL*(L1 + L2);
roL = mL/(L1 + L2);
Ei = 7.2*10^10;
Ai = bL*hL^3/12;
rouB = 2700;
lB = 0.15;
hB = 0.06;
tB = 0.08;
mB = rouB*lB*hB*tB;
JB = mB*(lB^2 + hB^2)/12;
k = 4.68*10^9;
mu = 0.48;
g = 9.8;
sita=3*pi/2;
syms q1t q2t q3t q5t q7t q1tp q2tp q3tp q5tp q7tp real
syms x1 x2 real
syms sigmu real
%defining functions
fe1=-(x1^2/L1^2) + x1^3/L1^3;
fec1=-((2*x1^3)/(3*L1^4));
fe2=-1/(15*L1) + x2/(15*L1*L2);
fe3=(-3*x2)/L2 + (6*x2^2)/L2^2 - (4*x2^3)/L2^3;
fec2q54x2=-x2/(450*L1^2*L2^2);
fec2q52q7x2=x2/(5*L1*L2^2) - (2*x2^2)/(5*L1*L2^3) + (4*x2^3)/(15*L1*L2^4);
fec2q72x2=(-9*x2)/(2*L2^2) + (18*x2^2)/L2^3 - (36*x2^3)/L2^4;
infe1x1=double(int(roL*fe1,x1,0,L1));
infe1x1pf=double(int(roL*fe1^2,x1,0,L1));
infe1x1fec1x1=double(int(roL*fe1*fec1,x1,0,L1));
inx1=double(int(roL*x1,x1,0,L1));
inx1pf=double(int(roL*x1^2,x1,0,L1));
inx1fe1x1=double(int(roL*x1*fe1,x1,0,L1));
inx1fec1x1=double(int(roL*x1*fec1,x1,0,L1));
infec1x1=double(int(roL*fec1,x1,0,L1));
infec1x1pf=double(int(roL*fec1^2,x1,0,L1));
intL1=roL*L1;
infe2x2=double(int(roL*fe2,x2,0,L2));
infe2x2pf=double(int(roL*fe2^2,x2,0,L2));
infe2x2fe3x2=double(int(roL*fe2*fe3,x2,0,L2));
infe2x2fec2q72x2=double(int(roL*fe2*fec2q72x2,x2,0,L2));
infe2x2fec2q52q7x2=double(int(roL*fe2*fec2q52q7x2,x2,0,L2));
infe2x2fec2q54x2=double(int(roL*fe2*fec2q54x2,x2,0,L2));
inx2fe2x2=double(int(roL*fe2*x2,x2,0,L2));
inx2=double(int(roL*x2,x2,0,L2));
inx2pf=double(int(roL*x2^2,x2,0,L2));
inx2fec2q52q7x2=double(int(roL*x2*fec2q52q7x2,x2,0,L2));
inx2fec2q54x2=double(int(roL*fec2q54x2*x2,x2,0,L2));
inx2fec2q72x2=double(int(roL*fec2q72x2*x2,x2,0,L2));
inx2fe3x2=double(int(roL*x2*fe3,x2,0,L2));
infe3x2=double(int(roL*fe3,x2,0,L2));
infe3x2pf=double(int(roL*fe3^2,x2,0,L2));
infe3x2fec2q52q7x2=double(int(roL*fe3*fec2q52q7x2,x2,0,L2));
infe3x2fec2q54x2=double(int(roL*fe3*fec2q54x2,x2,0,L2));
infe3x2fec2q72x2=double(int(roL*fe3*fec2q72x2,x2,0,L2));
infec2q52q7x2fec2q54x2=double(int(roL*fec2q52q7x2*fec2q54x2,x2,0,L2));
infec2q54x2=double(int(roL*fec2q54x2,x2,0,L2));
infec2q54x2fec2q72x2=double(int(roL*fec2q54x2*fec2q72x2,x2,0,L2));
infec2q54x2pf=double(int(roL*fec2q54x2^2,x2,0,L2));
infec2q72x2=double(int(roL*fec2q72x2,x2,0,L2));
infec2q72x2pf=double(int(roL*fec2q72x2^2,x2,0,L2));
infec2q52q7x2=double(int(roL*fec2q52q7x2,x2,0,L2));
infec2q52q7x2pf=double(int(roL*fec2q52q7x2^2,x2,0,L2));
infec2q52q7x2fec2q72x2=double(int(roL*fec2q52q7x2*fec2q72x2,x2,0,L2));
intL2=roL*L2;
q1tppf_TL1=intL1/2;
q2tppf_TL1=intL1/2;
q1tpq3tp_TL1=-(cos(q3t)*infe1x1*q5t + inx1*sin(q3t) + infec1x1*q5t^2*sin(q3t));
q2tpq3tp_TL1=inx1*cos(q3t) + cos(q3t)*infec1x1*q5t^2 - infe1x1*q5t*sin(q3t);
q3tppf_TL1=inx1pf/2 + infe1x1pf*q5t^2/2 + inx1fec1x1*q5t^2 + infec1x1pf*q5t^4/2;
q1tpq5tp_TL1=2*cos(q3t)*infec1x1*q5t - infe1x1*sin(q3t);
q2tpq5tp_TL1=cos(q3t)*infe1x1 + 2*infec1x1*q5t*sin(q3t);
q3tpq5tp_TL1=inx1fe1x1- infe1x1fec1x1*q5t^2;
q5tppf_TL1=infe1x1pf/2 + 2*infec1x1pf*q5t^2;
q1tppf_TL2=intL2/2;
q2tppf_TL2=intL2/2;
q1tpq3tp_TL2=infe3x2*q7t*sin(q3t)*sin(sita) - inx2*cos(sita)*sin(q3t) -...
L1*intL2*sin(q3t)-...
inx2*cos(q3t)*sin(sita) - infe2x2*q5t^2*cos(q3t)*cos(sita) - ...
infec2q54x2*q5t^4*cos(q3t)*sin(sita) - infec2q54x2*q5t^4*cos(sita)*sin(q3t) -...
infec2q72x2*q7t^2*cos(q3t)*sin(sita) - infec2q72x2*q7t^2*cos(sita)*sin(q3t) +...
infe2x2*q5t^2*sin(q3t)*sin(sita) - infe3x2*q7t*cos(q3t)*cos(sita) -...
infec2q52q7x2*q5t^2*q7t*cos(q3t)*sin(sita) - infec2q52q7x2*q5t^2*q7t*cos(sita)*sin(q3t);
q2tpq3tp_TL2=inx2*cos(q3t)*cos(sita) - inx2*sin(q3t)*sin(sita) + L1*intL2*cos(q3t) +...
infec2q54x2*q5t^4*cos(q3t)*cos(sita) + infec2q72x2*q7t^2*cos(q3t)*cos(sita) - ...
infe2x2*q5t^2*cos(q3t)*sin(sita) - infe2x2*q5t^2*cos(sita)*sin(q3t) -...
infec2q54x2*q5t^4*sin(q3t)*sin(sita) - infec2q72x2*q7t^2*sin(q3t)*sin(sita) -...
infe3x2*q7t*cos(q3t)*sin(sita) - infe3x2*q7t*cos(sita)*sin(q3t) + ...
infec2q52q7x2*q5t^2*q7t*cos(q3t)*cos(sita) - infec2q52q7x2*q5t^2*q7t*sin(q3t)*sin(sita);
q3tppf_TL2=inx2pf/2 + (L1^2*intL2)/2 + (infe2x2pf*q5t^4)/2 + (infec2q54x2pf*q5t^8)/2 +...
inx2fec2q54x2*q5t^4 + (infe3x2pf*q7t^2)/2 + (infec2q72x2pf*q7t^4)/2 +...
inx2fec2q72x2*q7t^2 + infe2x2fe3x2*q5t^2*q7t + infec2q52q7x2fec2q54x2*q5t^6*q7t +...
inx2fec2q52q7x2*q5t^2*q7t + ...
(infec2q52q7x2pf*q5t^4*q7t^2)/2 + infec2q54x2fec2q72x2*q5t^4*q7t^2 +...
infec2q52q7x2fec2q72x2*q5t^2*q7t^3 + L1*inx2*cos(sita) - L1*infe3x2*q7t*sin(sita) + ...
L1*infec2q54x2*q5t^4*cos(sita) + L1*infec2q72x2*q7t^2*cos(sita) -...
L1*infe2x2*q5t^2*sin(sita) + ...
L1*infec2q52q7x2*q5t^2*q7t*cos(sita);
q1tpq5tp_TL2=4*infec2q54x2*q5t^3*cos(q3t)*cos(sita) -...
4*infec2q54x2*q5t^3*sin(q3t)*sin(sita) -...
2*infe2x2*q5t*cos(q3t)*sin(sita) - 2*infe2x2*q5t*cos(sita)*sin(q3t) + ...
2*infec2q52q7x2*q5t*q7t*cos(q3t)*cos(sita) - 2*infec2q52q7x2*q5t*q7t*sin(q3t)*sin(sita);
q2tpq5tp_TL2=4*infec2q54x2*q5t^3*cos(q3t)*sin(sita) -...
2*infe2x2*q5t*sin(q3t)*sin(sita) +...
4*infec2q54x2*q5t^3*cos(sita)*sin(q3t) + 2*infe2x2*q5t*cos(q3t)*cos(sita) + ...
2*infec2q52q7x2*q5t*q7t*cos(q3t)*sin(sita) + 2*infec2q52q7x2*q5t*q7t*cos(sita)*sin(q3t);
q3tpq5tp_TL2=2*inx2fe2x2*q5t - 2*infe2x2fec2q54x2*q5t^5 +2*infe2x2fec2q72x2*q5t*q7t^2 ...
- 4*infe3x2fec2q54x2*q5t^3*q7t - 2*infe3x2fec2q52q7x2*q5t*q7t^2 + ...
2*L1*infe2x2*q5t*cos(sita) + 4*L1*infec2q54x2*q5t^3*sin(sita) +...
2*L1*infec2q52q7x2*q5t*q7t*sin(sita);
q5tppf_TL2=8*infec2q54x2pf*q5t^6 + 8*infec2q52q7x2fec2q54x2*q5t^4*q7t +...
2*infec2q52q7x2pf*q5t^2*q7t^2 +...
2*infe2x2pf*q5t^2;
q1tpq7tp_TL2=infec2q52q7x2*q5t^2*cos(q3t)*cos(sita) - infe3x2*cos(sita)*sin(q3t) -...
2*infec2q72x2*q7t*sin(q3t)*sin(sita) - infe3x2*cos(q3t)*sin(sita) - ...
infec2q52q7x2*q5t^2*sin(q3t)*sin(sita) + 2*infec2q72x2*q7t*cos(q3t)*cos(sita);
q2tpq7tp_TL2=infe3x2*cos(q3t)*cos(sita) - infe3x2*sin(q3t)*sin(sita) + ...
infec2q52q7x2*q5t^2*cos(q3t)*sin(sita) + infec2q52q7x2*q5t^2*cos(sita)*sin(q3t) +...
2*infec2q72x2*q7t*cos(q3t)*sin(sita) + 2*infec2q72x2*q7t*cos(sita)*sin(q3t);
q3tpq7tp_TL2=inx2fe3x2 + L1*cos(sita)*infe3x2 - infe2x2fec2q52q7x2*q5t^4 +...
infe3x2fec2q54x2*q5t^4 - 2*infe2x2fec2q72x2*q5t^2*q7t - infe3x2fec2q72x2*q7t^2 +...
L1*infec2q52q7x2*q5t^2*sin(sita) + 2*L1*infec2q72x2*q7t*sin(sita);
q5tpq7tp_TL2=2*infe2x2fe3x2*q5t + 4*infec2q52q7x2fec2q54x2*q5t^5 +...
2*infec2q52q7x2pf*q5t^3*q7t + ...
8*infec2q54x2fec2q72x2*q5t^3*q7t + 4*infec2q52q7x2fec2q72x2*q5t*q7t^2;
q7tppf_TL2=(infec2q52q7x2pf*q5t^4)/2 + 2*infec2q52q7x2fec2q72x2*q5t^2*q7t +...
2*infec2q72x2pf*q7t^2 + infe3x2pf/2;
q1tppf_TB=mB/2;
q2tppf_TB=mB/2;
q3tppf_TB=JB/2;
q1tppf=q1tppf_TL1+q1tppf_TL2+q1tppf_TB;
q2tppf=q2tppf_TL1+q2tppf_TL2+q2tppf_TB;
q1tpq3tp=q1tpq3tp_TL1+q1tpq3tp_TL2;
q2tpq3tp=q2tpq3tp_TL1+q2tpq3tp_TL2;
q3tppf=q3tppf_TL1+q3tppf_TL2+q3tppf_TB;
q1tpq5tp=q1tpq5tp_TL1+q1tpq5tp_TL2;
q2tpq5tp=q2tpq5tp_TL1+q2tpq5tp_TL2;
q3tpq5tp=q3tpq5tp_TL1+q3tpq5tp_TL2;
q5tppf=q5tppf_TL1+q5tppf_TL2;
q1tpq7tp=q1tpq7tp_TL2;
q2tpq7tp=q2tpq7tp_TL2;
q3tpq7tp=q3tpq7tp_TL2;
q5tpq7tp=q5tpq7tp_TL2;
q7tppf=q7tppf_TL2;
M11=2*q1tppf;
M12=0;
M13=q1tpq3tp;
M15=q1tpq5tp;
M17=q1tpq7tp;
M21=M12;
M22=2*q2tppf;
M23=q2tpq3tp;
M25=q2tpq5tp;
M27=q2tpq7tp;
M31=M13;
M32=M23;
M33=2*q3tppf;
M35=q3tpq5tp;
M37=q3tpq7tp;
M51=M15;
M52=M25;
M53=M35;
M55=2*q5tppf;
M57=q5tpq7tp;
M71=M11;
M72=M27;
M73=M37;
M75=M57;
M77=2*q7tppf;
fenmu=(M77*M12^2*M35^2 - 2*M12^2*M35*M37*M57 + M55*M12^2*M37^2 + M33*M12^2*M57^2 - ...
M33*M55*M77*M12^2 - 2*M12*M13*M23*M57^2 + 2*M55*M77*M12*M13*M23 -...
2*M77*M12*M13*M25*M35 + 2*M12*M13*M25*M37*M57 + 2*M12*M13*M27*M35*M57 - ...
2*M55*M12*M13*M27*M37 - 2*M77*M12*M15*M23*M35 + 2*M12*M15*M23*M37*M57 -...
2*M12*M15*M25*M37^2 + 2*M33*M77*M12*M15*M25 + 2*M12*M15*M27*M35*M37 - ...
2*M33*M12*M15*M27*M57 + 2*M12*M17*M23*M35*M57 - 2*M55*M12*M17*M23*M37 +...
2*M12*M17*M25*M35*M37 - 2*M33*M12*M17*M25*M57 - 2*M12*M17*M27*M35^2 + ...
2*M33*M55*M12*M17*M27 + M77*M13^2*M25^2 - 2*M13^2*M25*M27*M57 +...
M55*M13^2*M27^2 + M22*M13^2*M57^2 - M22*M55*M77*M13^2 - 2*M77*M13*M15*M23*M25 +...
2*M13*M15*M23*M27*M57 + 2*M13*M15*M25*M27*M37 - 2*M13*M15*M27^2*M35 + ...
2*M22*M77*M13*M15*M35 - 2*M22*M13*M15*M37*M57 + 2*M13*M17*M23*M25*M57 - ...
2*M55*M13*M17*M23*M27 - 2*M13*M17*M25^2*M37 + 2*M13*M17*M25*M27*M35 - ...
2*M22*M13*M17*M35*M57 + 2*M22*M55*M13*M17*M37 + M77*M15^2*M23^2 - ...
2*M15^2*M23*M27*M37 + M33*M15^2*M27^2 + M22*M15^2*M37^2 - M22*M33*M77*M15^2 - ...
2*M15*M17*M23^2*M57 + 2*M15*M17*M23*M25*M37 + 2*M15*M17*M23*M27*M35 -...
2*M33*M15*M17*M25*M27 - 2*M22*M15*M17*M35*M37 + 2*M22*M33*M15*M17*M57 +...
M55*M17^2*M23^2 - 2*M17^2*M23*M25*M35 + M33*M17^2*M25^2 + M22*M17^2*M35^2 -...
M22*M33*M55*M17^2 + M11*M23^2*M57^2 - M11*M55*M77*M23^2 + 2*M11*M77*M23*M25*M35 - ...
2*M11*M23*M25*M37*M57 - 2*M11*M23*M27*M35*M57 + 2*M11*M55*M23*M27*M37 + ...
M11*M25^2*M37^2 - M11*M33*M77*M25^2 - 2*M11*M25*M27*M35*M37 + 2*M11*M33*M25*M27*M57 +...
M11*M27^2*M35^2 - M11*M33*M55*M27^2 - M11*M22*M77*M35^2 +2*M11*M22*M35*M37*M57-...
M11*M22*M55*M37^2 -...
M11*M22*M33*M57^2 + M11*M22*M33*M55*M77);
InvM11=-(- M23^2*M57^2 + M55*M77*M23^2 - 2*M77*M23*M25*M35 + 2*M23*M25*M37*M57 +...
2*M23*M27*M35*M57 - 2*M55*M23*M27*M37 - M25^2*M37^2 + M33*M77*M25^2 +...
2*M25*M27*M35*M37 - 2*M33*M25*M27*M57 - M27^2*M35^2 + M33*M55*M27^2 +...
M22*M77*M35^2 - 2*M22*M35*M37*M57 +...
M22*M55*M37^2 + M22*M33*M57^2 - M22*M33*M55*M77)/fenmu;
InvM12=(M12*M33*M57^2 - M17*M27*M35^2 - M13*M23*M57^2 - M15*M25*M37^2 + M12*M37^2*M55+...
M12*M35^2*M77 + M15*M27*M35*M37 + M17*M25*M35*M37 + M13*M25*M37*M57 + M13*M27*M35*M57-...
M13*M27*M37*M55 + M15*M23*M37*M57 - M15*M27*M33*M57 + M17*M23*M35*M57 -...
M17*M23*M37*M55 - M17*M25*M33*M57 + M17*M27*M33*M55 - 2*M12*M35*M37*M57 -...
M13*M25*M35*M77 - M15*M23*M35*M77+ M15*M25*M33*M77 +...
M13*M23*M55*M77 -M12*M33*M55*M77)/fenmu;
InvM13=(M13*M22*M57^2 - M17*M25^2*M37 - M12*M23*M57^2 - M15*M27^2*M35 + M13*M27^2*M55+...
M13*M25^2*M77 + M15*M25*M27*M37 + M17*M25*M27*M35 - 2*M13*M25*M27*M57 +...
M15*M23*M27*M57 +M17*M23*M25*M57 - M17*M23*M27*M55 + M12*M25*M37*M57 +...
M12*M27*M35*M57 - M12*M27*M37*M55 -M15*M22*M37*M57 - M17*M22*M35*M57 +...
M17*M22*M37*M55 - M15*M23*M25*M77 - M12*M25*M35*M77 +...
M15*M22*M35*M77 + M12*M23*M55*M77 - M13*M22*M55*M77)/fenmu;
InvM15=(M15*M22*M37^2 - M12*M25*M37^2 - M13*M27^2*M35 + M15*M27^2*M33 -M17*M23^2*M57 +...
M15*M23^2*M77 +M13*M25*M27*M37 -2*M15*M23*M27*M37 +M17*M23*M25*M37 + M17*M23*M27*M35 -...
M17*M25*M27*M33 + M12*M27*M35*M37 - M17*M22*M35*M37 +M13*M23*M27*M57 +M12*M23*M37*M57-...
M12*M27*M33*M57-M13*M22*M37*M57+M17*M22*M33*M57-M13*M23*M25*M77 - M12*M23*M35*M77 +...
M12*M25*M33*M77 + M13*M22*M35*M77 - M15*M22*M33*M77)/fenmu;
InvM17=(M17*M22*M35^2 - M12*M27*M35^2 - M13*M25^2*M37 + M17*M25^2*M33 - M15*M23^2*M57+...
M17*M23^2*M55 + M13*M25*M27*M35 + M15*M23*M25*M37 + M15*M23*M27*M35 - M15*M25*M27*M33-...
2*M17*M23*M25*M35 + M12*M25*M35*M37 -M15*M22*M35*M37 +M13*M23*M25*M57 -...
M13*M23*M27*M55+ M12*M23*M35*M57 - M12*M23*M37*M55 - M12*M25*M33*M57 +...
M12*M27*M33*M55 - M13*M22*M35*M57 + M13*M22*M37*M55 +...
M15*M22*M33*M57 - M17*M22*M33*M55)/fenmu;
InvM21=InvM12;
InvM22=-(- M13^2*M57^2 + M55*M77*M13^2 - 2*M77*M13*M15*M35 + 2*M13*M15*M37*M57 +...
2*M13*M17*M35*M57 - 2*M55*M13*M17*M37 - M15^2*M37^2 + M33*M77*M15^2 +...
2*M15*M17*M35*M37 -2*M33*M15*M17*M57 - M17^2*M35^2 + M33*M55*M17^2 +...
M11*M77*M35^2 - 2*M11*M35*M37*M57 + M11*M55*M37^2 + M11*M33*M57^2 -...
M11*M33*M55*M77)/fenmu;
InvM23=(M11*M23*M57^2 - M15^2*M27*M37 - M12*M13*M57^2 - M17^2*M25*M35 +...
M17^2*M23*M55 + M15^2*M23*M77 + M15*M17*M25*M37 + M15*M17*M27*M35 + ...
M13*M15*M27*M57 + M13*M17*M25*M57 - M13*M17*M27*M55 - 2*M15*M17*M23*M57 +...
M12*M15*M37*M57 + M12*M17*M35*M57 - M12*M17*M37*M55 - M11*M25*M37*M57 -...
M11*M27*M35*M57 + M11*M27*M37*M55 - M13*M15*M25*M77 - M12*M15*M35*M77+ ...
M11*M25*M35*M77 + M12*M13*M55*M77 - M11*M23*M55*M77)/fenmu;
InvM25=(M11*M25*M37^2 - M12*M15*M37^2 - M17^2*M23*M35 + M17^2*M25*M33 -...
M13^2*M27*M57 + M13^2*M25*M77 + M13*M15*M27*M37 - 2*M13*M17*M25*M37 +...
M13*M17*M27*M35 + M15*M17*M23*M37 - M15*M17*M27*M33 + M12*M17*M35*M37 -...
M11*M27*M35*M37 + M13*M17*M23*M57 + M12*M13*M37*M57 -...
M12*M17*M33*M57 - M11*M23*M37*M57 + M11*M27*M33*M57 - M13*M15*M23*M77 -...
M12*M13*M35*M77 +M12*M15*M33*M77 + M11*M23*M35*M77 - M11*M25*M33*M77)/fenmu;
InvM27=(M11*M27*M35^2 - M12*M17*M35^2 - M15^2*M23*M37 + M15^2*M27*M33 -...
M13^2*M25*M57 + M13^2*M27*M55 + M13*M15*M25*M37 - 2*M13*M15*M27*M35 +...
M13*M17*M25*M35 + M15*M17*M23*M35 - M15*M17*M25*M33 + M12*M15*M35*M37 -...
M11*M25*M35*M37 + M13*M15*M23*M57 - M13*M17*M23*M55 +...
M12*M13*M35*M57 - M12*M13*M37*M55 - M12*M15*M33*M57 + M12*M17*M33*M55 -...
M11*M23*M35*M57 +M11*M23*M37*M55 + M11*M25*M33*M57 - M11*M27*M33*M55)/fenmu;
InvM31=InvM13;
InvM32=InvM23;
InvM33=-(- M12^2*M57^2 + M55*M77*M12^2 - 2*M77*M12*M15*M25 + 2*M12*M15*M27*M57 +...
2*M12*M17*M25*M57 - 2*M55*M12*M17*M27 - M15^2*M27^2 + M22*M77*M15^2 +...
2*M15*M17*M25*M27 - 2*M22*M15*M17*M57 - M17^2*M25^2 + M22*M55*M17^2 +...
M11*M77*M25^2 - 2*M11*M25*M27*M57 + M11*M55*M27^2 + M11*M22*M57^2 -...
M11*M22*M55*M77)/fenmu;
InvM35=(M11*M27^2*M35 - M17^2*M23*M25 - M13*M15*M27^2 + M17^2*M22*M35 -...
M12^2*M37*M57 +M12^2*M35*M77 + M13*M17*M25*M27 + M15*M17*M23*M27 + ...
M12*M15*M27*M37 + M12*M17*M25*M37 -2*M12*M17*M27*M35 -M15*M17*M22*M37- ...
M11*M25*M27*M37 + M12*M13*M27*M57 + M12*M17*M23*M57- M13*M17*M22*M57 -...
M11*M23*M27*M57 + M11*M22*M37*M57 - M12*M13*M25*M77 - M12*M15*M23*M77+...
M13*M15*M22*M77 + M11*M23*M25*M77 - M11*M22*M35*M77)/fenmu;
InvM37=(M11*M25^2*M37 - M15^2*M23*M27 - M13*M17*M25^2 + M15^2*M22*M37 -M12^2*M35*M57 +...
M12^2*M37*M55 + M13*M15*M25*M27 + M15*M17*M23*M25 - 2*M12*M15*M25*M37 +...
M12*M15*M27*M35 + M12*M17*M25*M35 - M15*M17*M22*M35 - M11*M25*M27*M35 + ...
M12*M13*M25*M57 - M12*M13*M27*M55 + M12*M15*M23*M57 - M12*M17*M23*M55 -...
M13*M15*M22*M57 + M13*M17*M22*M55 - M11*M23*M25*M57 +...
M11*M23*M27*M55 + M11*M22*M35*M57 - M11*M22*M37*M55)/fenmu;
InvM51=InvM15;
InvM52=InvM25;
InvM53=InvM35;
InvM55=-(- M12^2*M37^2 + M33*M77*M12^2 - 2*M77*M12*M13*M23 + 2*M12*M13*M27*M37 +...
2*M12*M17*M23*M37 - 2*M33*M12*M17*M27 - M13^2*M27^2 + M22*M77*M13^2+ ...
2*M13*M17*M23*M27 - 2*M22*M13*M17*M37 - M17^2*M23^2 + M22*M33*M17^2 + ...
M11*M77*M23^2 - 2*M11*M23*M27*M37 + M11*M33*M27^2 + M11*M22*M37^2 -...
M11*M22*M33*M77)/fenmu;
InvM57=(M11*M23^2*M57 - M13^2*M25*M27 - M12^2*M35*M37 - M15*M17*M23^2 + ...
M13^2*M22*M57 +M12^2*M33*M57 + M13*M15*M23*M27 + M13*M17*M23*M25 +...
M12*M13*M25*M37 + M12*M13*M27*M35 + M12*M15*M23*M37 - M12*M15*M27*M33 +...
M12*M17*M23*M35 - M12*M17*M25*M33 - M13*M15*M22*M37 - M13*M17*M22*M35 +...
M15*M17*M22*M33 - M11*M23*M25*M37 - M11*M23*M27*M35 + M11*M25*M27*M33 +...
M11*M22*M35*M37 - 2*M12*M13*M23*M57 - M11*M22*M33*M57)/fenmu;
InvM71=InvM17;
InvM72=InvM27;
InvM73=InvM37;
InvM75=InvM57;
InvM77=-(- M12^2*M35^2 + M33*M55*M12^2 - 2*M55*M12*M13*M23 + 2*M12*M13*M25*M35 +...
2*M12*M15*M23*M35 - 2*M33*M12*M15*M25 - M13^2*M25^2 + M22*M55*M13^2 + ...
2*M13*M15*M23*M25 - 2*M22*M13*M15*M35 - M15^2*M23^2 + M22*M33*M15^2 +...
M11*M55*M23^2 - 2*M11*M23*M25*M35 + M11*M33*M25^2 + M11*M22*M35^2 -...
M11*M22*M33*M55)/fenmu;
InvMty22=[InvM11,InvM12,InvM13,InvM15,InvM17;InvM21,InvM22,InvM23,InvM25,InvM27;...
InvM31,InvM32,InvM33,InvM35,InvM37;InvM51,InvM52,InvM53,InvM55,InvM57;...
InvM71,InvM72,InvM73,InvM75,InvM77];
%elastic energy parameters
fe1pp=diff(fe1,x1,2);
intkff1=fe1pp^2;
kff1=double(int(Ai*Ei*intkff1,x1,0,L1));
fe2pp=diff(fe2,x2,2);
intkff20=fe2pp^2;
kff2_2pf=double(int(Ai*Ei*intkff20,x2,0,L2));
fe3pp=diff(fe3,x2,2);
intkff22=fe3pp^2;
kff2_3pf=double(int(Ai*Ei*intkff22,x2,0,L2));
intkff23=fe2pp*fe3pp;
kff2_23=double(int(Ai*Ei*intkff23,x2,0,L2));
%elastic energy
PeL1=Ai*Ei*q5t^2*kff1/2;
PeL2=Ai*Ei*(q5t^4*kff2_2pf+2*q5t^2*q7t*kff2_23+q7t^2*kff2_3pf)/2;
PeL=PeL1+PeL2;
%gravitational energy of L
PLg1=(q2t*roL*L1 + cos (q3t)*infe1x1*q5t + inx1*sin (q3t) +...
infec1x1*q5t^2*sin (q3t))*g;
PLg2=g*(q2t*L2*roL + cos (sita)*cos (q3t)*infe2x2*q5t^2 + ...
cos (sita)*cos (q3t)*infe3x2*q7t + inx2*cos (q3t)*sin (sita) + ...
cos (q3t)*infec2q54x2*q5t^4*sin (sita) + ...
cos (q3t)*infec2q52q7x2*q5t^2*q7t*sin (sita) + ...
cos (q3t)*infec2q72x2*q7t^2*sin (sita) + L1*sin (q3t)*L2*roL +...
inx2*cos (sita)*sin (q3t) + cos (sita)*infec2q54x2*q5t^4*sin (q3t) + ...
cos (sita)*infec2q52q7x2*q5t^2*q7t*sin (q3t) + ...
cos (sita)*infec2q72x2*q7t^2*sin (q3t) - ...
infe2x2*q5t^2*sin (sita)*sin (q3t) - infe3x2*q7t*sin (sita)*sin (q3t));
PLg=PLg1+PLg2;
%gravitational energy of B
PBg=mB*g*q2t;
%gravitational energy of R
PRg=PLg+PBg;
fec2q54L2=subs(fec2q54x2,x2,L2);
fec2q52q7L2=subs(fec2q52q7x2,x2,L2);
fec2q72L2=subs(fec2q72x2,x2,L2);
fe2L2=subs(fe2,x2,L2);
fe3L2=subs(fe3,x2,L2);
r0pEx=L1*cos (q3t) + L2*cos (sita)*cos (q3t) + q1t + ...
cos (sita)*cos (q3t)*fec2q54L2*q5t^4 + ...
cos (sita)*cos (q3t)*fec2q52q7L2*q5t^2*q7t + ...
cos (sita)*cos (q3t)*fec2q72L2*q7t^2 - ...
cos (q3t)*fe2L2*q5t^2*sin (sita) - cos (q3t)*fe3L2*q7t*sin (sita) - ...
cos (sita)*fe2L2*q5t^2*sin (q3t) - cos (sita)*fe3L2*q7t*sin (q3t) -...
L2*sin (sita)*sin (q3t) - fec2q54L2*q5t^4*sin (sita)*sin (q3t) - ...
fec2q52q7L2*q5t^2*q7t*sin (sita)*sin (q3t) - ...
fec2q72L2*q7t^2*sin (sita)*sin (q3t);
r0pEy=q2t + cos (sita)*cos (q3t)*fe2L2*q5t^2 +...
cos (sita)*cos (q3t)*fe3L2*q7t + L2*cos (q3t)*sin (sita) +...
cos (q3t)*fec2q54L2*q5t^4*sin (sita) +...
cos (q3t)*fec2q52q7L2*q5t^2*q7t*sin (sita) + ...
cos (q3t)*fec2q72L2*q7t^2*sin (sita) + L1*sin (q3t) +...
L2*cos (sita)*sin (q3t) + cos (sita)*fec2q54L2*q5t^4*sin (q3t) + ...
cos (sita)*fec2q52q7L2*q5t^2*q7t*sin (q3t) + ...
cos (sita)*fec2q72L2*q7t^2*sin (q3t) - ...
fe2L2*q5t^2*sin (sita)*sin (q3t) - fe3L2*q7t*sin (sita)*sin (q3t);
r0pE=[r0pEx;r0pEy];
vpEx=q1tp - cos (sita)*cos (q3t)*fe2L2*q5t^2*q3tp -...
cos (sita)*cos (q3t)*fe3L2*q7t*q3tp - L2*cos (q3t)*sin (sita)*q3tp -...
cos (q3t)*fec2q54L2*q5t^4*sin (sita)*q3tp - ...
cos (q3t)*fec2q52q7L2*q5t^2*q7t*sin (sita)*q3tp - ...
cos (q3t)*fec2q72L2*q7t^2*sin (sita)*q3tp - L1*sin (q3t)*q3tp -...
L2*cos (sita)*sin (q3t)*q3tp - ...
cos (sita)*fec2q54L2*q5t^4*sin (q3t)*q3tp - ...
cos (sita)*fec2q52q7L2*q5t^2*q7t*sin (q3t)*q3tp - ...
cos (sita)*fec2q72L2*q7t^2*sin (q3t)*q3tp + ...
fe2L2*q5t^2*sin (sita)*sin (q3t)*q3tp + ...
fe3L2*q7t*sin (sita)*sin (q3t)*q3tp + ...
4*cos (sita)*cos (q3t)*fec2q54L2*q5t^3*q5tp + ...
2*cos (sita)*cos (q3t)*fec2q52q7L2*q5t*q7t*q5tp - ...
2*cos (q3t)*fe2L2*q5t*sin (sita)*q5tp - ...
2*cos (sita)*fe2L2*q5t*sin (q3t)*q5tp - ...
4*fec2q54L2*q5t^3*sin (sita)*sin (q3t)*q5tp - ...
2*fec2q52q7L2*q5t*q7t*sin (sita)*sin (q3t)*q5tp + ...
cos (sita)*cos (q3t)*fec2q52q7L2*q5t^2*q7tp + ...
2*cos (sita)*cos (q3t)*fec2q72L2*q7t*q7tp - ...
cos (q3t)*fe3L2*sin (sita)*q7tp - cos (sita)*fe3L2*sin (q3t)*q7tp - ...
fec2q52q7L2*q5t^2*sin (sita)*sin (q3t)*q7tp - ...
2*fec2q72L2*q7t*sin (sita)*sin (q3t)*q7tp;
RLEq1=diff(r0pE,q1t);
RLEq2=diff(r0pE,q2t);
RLEq3=diff(r0pE,q3t);
RLEq5=diff(r0pE,q5t);
RLEq7=diff(r0pE,q7t);
Fn=k*abs(r0pEy)^(3/2);
Ft=sigmu*mu*Fn;
F=[Ft;Fn];
Q1=RLEq1.'*F;
Q2=RLEq2.'*F;
Q3=RLEq3.'*F;
Q5=RLEq5.'*F;
Q7=RLEq7.'*F;
Q=[Q1;Q2;Q3;Q5;Q7];
infe1x1pppf=kff1;
infe2x2pppf=kff2_2pf;
infe2x2ppfe3x2pp=kff2_23;
infe3x2pppf=kff2_3pf;
the left part of this program is posted here enter link description here
