How to accelerate the following code using 'equationsToMatrix'?

i just want to solve the 20 linear equations with 'equationsToMatrix',but it really runs out a lot of time,so how to accelerate it?

codes are as this

syms a0 a1 a2 a3 a4  a5 a6 a7 a8 a9 alpha0 alpha1 alpha2 alpha3 alpha4 alpha5 alpha6 alpha7 alpha8 alpha9; 
 
omega=2.000000;
  
m0=0.685725;
m1=3.006665;
m2=6.217633;
m3=9.381297;
m4=12.533817;
m5=15.681942;
m6=18.827881;
m7=21.972575;
m8=25.116492;
m9=28.259892;
m10=31.402930;

h=1;
 
hgang=1;
 
agang=5;
 
dgang=1;
 
 
a=agang/hgang;
 
d=dgang/hgang; 
 
equation0=...
 ...
alpha0 == (2*a1*sin(d*m1)*besselk(1, a*m1))/(d*m1*(sin(2*m1)/(4*m1) + 1/2)^(1/2)) + (2*a2*sin(d*m2)*besselk(1, a*m2))/(d*m2*(sin(2*m2)/(4*m2) + 1/2)^(1/2)) + (2*a3*sin(d*m3)*besselk(1, a*m3))/(d*m3*(sin(2*m3)/(4*m3) + 1/2)^(1/2)) + (2*a4*sin(d*m4)*besselk(1, a*m4))/(d*m4*(sin(2*m4)/(4*m4) + 1/2)^(1/2)) + (2*a5*sin(d*m5)*besselk(1, a*m5))/(d*m5*(sin(2*m5)/(4*m5) + 1/2)^(1/2)) + (2*a6*sin(d*m6)*besselk(1, a*m6))/(d*m6*(sin(2*m6)/(4*m6) + 1/2)^(1/2)) + (2*a7*sin(d*m7)*besselk(1, a*m7))/(d*m7*(sin(2*m7)/(4*m7) + 1/2)^(1/2)) + (2*a8*sin(d*m8)*besselk(1, a*m8))/(d*m8*(sin(2*m8)/(4*m8) + 1/2)^(1/2)) + (2*a9*sin(d*m9)*besselk(1, a*m9))/(d*m9*(sin(2*m9)/(4*m9) + 1/2)^(1/2)) + (2*a0*sinh(d*m0)*(besselj(1, a*m0) + bessely(1, a*m0)*i))/(d*m0*(sinh(2*m0)/(4*m0) + 1/2)^(1/2))
 
 
equation1=...
 ...
alpha1 == (2*a1*d*m1*sin(d*m1)*besselk(1, a*m1))/((sin(2*m1)/(4*m1) + 1/2)^(1/2)*(- d^2*m1^2 + pi^2)) + (2*a2*d*m2*sin(d*m2)*besselk(1, a*m2))/((sin(2*m2)/(4*m2) + 1/2)^(1/2)*(- d^2*m2^2 + pi^2)) + (2*a3*d*m3*sin(d*m3)*besselk(1, a*m3))/((sin(2*m3)/(4*m3) + 1/2)^(1/2)*(- d^2*m3^2 + pi^2)) + (2*a4*d*m4*sin(d*m4)*besselk(1, a*m4))/((sin(2*m4)/(4*m4) + 1/2)^(1/2)*(- d^2*m4^2 + pi^2)) + (2*a5*d*m5*sin(d*m5)*besselk(1, a*m5))/((sin(2*m5)/(4*m5) + 1/2)^(1/2)*(- d^2*m5^2 + pi^2)) + (2*a6*d*m6*sin(d*m6)*besselk(1, a*m6))/((sin(2*m6)/(4*m6) + 1/2)^(1/2)*(- d^2*m6^2 + pi^2)) + (2*a7*d*m7*sin(d*m7)*besselk(1, a*m7))/((sin(2*m7)/(4*m7) + 1/2)^(1/2)*(- d^2*m7^2 + pi^2)) + (2*a8*d*m8*sin(d*m8)*besselk(1, a*m8))/((sin(2*m8)/(4*m8) + 1/2)^(1/2)*(- d^2*m8^2 + pi^2)) + (2*a9*d*m9*sin(d*m9)*besselk(1, a*m9))/((sin(2*m9)/(4*m9) + 1/2)^(1/2)*(- d^2*m9^2 + pi^2)) - (2*a0*d*m0*sinh(d*m0)*(besselj(1, a*m0) + bessely(1, a*m0)*i))/((sinh(2*m0)/(4*m0) + 1/2)^(1/2)*(d^2*m0^2 + pi^2))
 
 
equation2=...
 ...
alpha2 == - (2*a1*d*m1*sin(d*m1)*besselk(1, a*m1))/((- d^2*m1^2 + 4*pi^2)*(sin(2*m1)/(4*m1) + 1/2)^(1/2)) - (2*a2*d*m2*sin(d*m2)*besselk(1, a*m2))/((- d^2*m2^2 + 4*pi^2)*(sin(2*m2)/(4*m2) + 1/2)^(1/2)) - (2*a3*d*m3*sin(d*m3)*besselk(1, a*m3))/((- d^2*m3^2 + 4*pi^2)*(sin(2*m3)/(4*m3) + 1/2)^(1/2)) - (2*a4*d*m4*sin(d*m4)*besselk(1, a*m4))/((- d^2*m4^2 + 4*pi^2)*(sin(2*m4)/(4*m4) + 1/2)^(1/2)) - (2*a5*d*m5*sin(d*m5)*besselk(1, a*m5))/((- d^2*m5^2 + 4*pi^2)*(sin(2*m5)/(4*m5) + 1/2)^(1/2)) - (2*a6*d*m6*sin(d*m6)*besselk(1, a*m6))/((- d^2*m6^2 + 4*pi^2)*(sin(2*m6)/(4*m6) + 1/2)^(1/2)) - (2*a7*d*m7*sin(d*m7)*besselk(1, a*m7))/((- d^2*m7^2 + 4*pi^2)*(sin(2*m7)/(4*m7) + 1/2)^(1/2)) - (2*a8*d*m8*sin(d*m8)*besselk(1, a*m8))/((- d^2*m8^2 + 4*pi^2)*(sin(2*m8)/(4*m8) + 1/2)^(1/2)) - (2*a9*d*m9*sin(d*m9)*besselk(1, a*m9))/((- d^2*m9^2 + 4*pi^2)*(sin(2*m9)/(4*m9) + 1/2)^(1/2)) + (2*a0*d*m0*sinh(d*m0)*(besselj(1, a*m0) + bessely(1, a*m0)*i))/((d^2*m0^2 + 4*pi^2)*(sinh(2*m0)/(4*m0) + 1/2)^(1/2))
 
 
equation3=...
 ...
alpha3 == (2*a1*d*m1*sin(d*m1)*besselk(1, a*m1))/((- d^2*m1^2 + 9*pi^2)*(sin(2*m1)/(4*m1) + 1/2)^(1/2)) + (2*a2*d*m2*sin(d*m2)*besselk(1, a*m2))/((- d^2*m2^2 + 9*pi^2)*(sin(2*m2)/(4*m2) + 1/2)^(1/2)) + (2*a3*d*m3*sin(d*m3)*besselk(1, a*m3))/((- d^2*m3^2 + 9*pi^2)*(sin(2*m3)/(4*m3) + 1/2)^(1/2)) + (2*a4*d*m4*sin(d*m4)*besselk(1, a*m4))/((- d^2*m4^2 + 9*pi^2)*(sin(2*m4)/(4*m4) + 1/2)^(1/2)) + (2*a5*d*m5*sin(d*m5)*besselk(1, a*m5))/((- d^2*m5^2 + 9*pi^2)*(sin(2*m5)/(4*m5) + 1/2)^(1/2)) + (2*a6*d*m6*sin(d*m6)*besselk(1, a*m6))/((- d^2*m6^2 + 9*pi^2)*(sin(2*m6)/(4*m6) + 1/2)^(1/2)) + (2*a7*d*m7*sin(d*m7)*besselk(1, a*m7))/((- d^2*m7^2 + 9*pi^2)*(sin(2*m7)/(4*m7) + 1/2)^(1/2)) + (2*a8*d*m8*sin(d*m8)*besselk(1, a*m8))/((- d^2*m8^2 + 9*pi^2)*(sin(2*m8)/(4*m8) + 1/2)^(1/2)) + (2*a9*d*m9*sin(d*m9)*besselk(1, a*m9))/((- d^2*m9^2 + 9*pi^2)*(sin(2*m9)/(4*m9) + 1/2)^(1/2)) - (2*a0*d*m0*sinh(d*m0)*(besselj(1, a*m0) + bessely(1, a*m0)*i))/((d^2*m0^2 + 9*pi^2)*(sinh(2*m0)/(4*m0) + 1/2)^(1/2))
 
 
equation4=...
 ...
alpha4 == - (2*a1*d*m1*sin(d*m1)*besselk(1, a*m1))/((- d^2*m1^2 + 16*pi^2)*(sin(2*m1)/(4*m1) + 1/2)^(1/2)) - (2*a2*d*m2*sin(d*m2)*besselk(1, a*m2))/((- d^2*m2^2 + 16*pi^2)*(sin(2*m2)/(4*m2) + 1/2)^(1/2)) - (2*a3*d*m3*sin(d*m3)*besselk(1, a*m3))/((- d^2*m3^2 + 16*pi^2)*(sin(2*m3)/(4*m3) + 1/2)^(1/2)) - (2*a4*d*m4*sin(d*m4)*besselk(1, a*m4))/((- d^2*m4^2 + 16*pi^2)*(sin(2*m4)/(4*m4) + 1/2)^(1/2)) - (2*a5*d*m5*sin(d*m5)*besselk(1, a*m5))/((- d^2*m5^2 + 16*pi^2)*(sin(2*m5)/(4*m5) + 1/2)^(1/2)) - (2*a6*d*m6*sin(d*m6)*besselk(1, a*m6))/((- d^2*m6^2 + 16*pi^2)*(sin(2*m6)/(4*m6) + 1/2)^(1/2)) - (2*a7*d*m7*sin(d*m7)*besselk(1, a*m7))/((- d^2*m7^2 + 16*pi^2)*(sin(2*m7)/(4*m7) + 1/2)^(1/2)) - (2*a8*d*m8*sin(d*m8)*besselk(1, a*m8))/((- d^2*m8^2 + 16*pi^2)*(sin(2*m8)/(4*m8) + 1/2)^(1/2)) - (2*a9*d*m9*sin(d*m9)*besselk(1, a*m9))/((- d^2*m9^2 + 16*pi^2)*(sin(2*m9)/(4*m9) + 1/2)^(1/2)) + (2*a0*d*m0*sinh(d*m0)*(besselj(1, a*m0) + bessely(1, a*m0)*i))/((d^2*m0^2 + 16*pi^2)*(sinh(2*m0)/(4*m0) + 1/2)^(1/2))
 
 
equation5=...
 ...
alpha5 == (2*a1*d*m1*sin(d*m1)*besselk(1, a*m1))/((- d^2*m1^2 + 25*pi^2)*(sin(2*m1)/(4*m1) + 1/2)^(1/2)) + (2*a2*d*m2*sin(d*m2)*besselk(1, a*m2))/((- d^2*m2^2 + 25*pi^2)*(sin(2*m2)/(4*m2) + 1/2)^(1/2)) + (2*a3*d*m3*sin(d*m3)*besselk(1, a*m3))/((- d^2*m3^2 + 25*pi^2)*(sin(2*m3)/(4*m3) + 1/2)^(1/2)) + (2*a4*d*m4*sin(d*m4)*besselk(1, a*m4))/((- d^2*m4^2 + 25*pi^2)*(sin(2*m4)/(4*m4) + 1/2)^(1/2)) + (2*a5*d*m5*sin(d*m5)*besselk(1, a*m5))/((- d^2*m5^2 + 25*pi^2)*(sin(2*m5)/(4*m5) + 1/2)^(1/2)) + (2*a6*d*m6*sin(d*m6)*besselk(1, a*m6))/((- d^2*m6^2 + 25*pi^2)*(sin(2*m6)/(4*m6) + 1/2)^(1/2)) + (2*a7*d*m7*sin(d*m7)*besselk(1, a*m7))/((- d^2*m7^2 + 25*pi^2)*(sin(2*m7)/(4*m7) + 1/2)^(1/2)) + (2*a8*d*m8*sin(d*m8)*besselk(1, a*m8))/((- d^2*m8^2 + 25*pi^2)*(sin(2*m8)/(4*m8) + 1/2)^(1/2)) + (2*a9*d*m9*sin(d*m9)*besselk(1, a*m9))/((- d^2*m9^2 + 25*pi^2)*(sin(2*m9)/(4*m9) + 1/2)^(1/2)) - (2*a0*d*m0*sinh(d*m0)*(besselj(1, a*m0) + bessely(1, a*m0)*i))/((d^2*m0^2 + 25*pi^2)*(sinh(2*m0)/(4*m0) + 1/2)^(1/2))
 
 
equation6=...
 ...
alpha6 == - (2*a1*d*m1*sin(d*m1)*besselk(1, a*m1))/((- d^2*m1^2 + 36*pi^2)*(sin(2*m1)/(4*m1) + 1/2)^(1/2)) - (2*a2*d*m2*sin(d*m2)*besselk(1, a*m2))/((- d^2*m2^2 + 36*pi^2)*(sin(2*m2)/(4*m2) + 1/2)^(1/2)) - (2*a3*d*m3*sin(d*m3)*besselk(1, a*m3))/((- d^2*m3^2 + 36*pi^2)*(sin(2*m3)/(4*m3) + 1/2)^(1/2)) - (2*a4*d*m4*sin(d*m4)*besselk(1, a*m4))/((- d^2*m4^2 + 36*pi^2)*(sin(2*m4)/(4*m4) + 1/2)^(1/2)) - (2*a5*d*m5*sin(d*m5)*besselk(1, a*m5))/((- d^2*m5^2 + 36*pi^2)*(sin(2*m5)/(4*m5) + 1/2)^(1/2)) - (2*a6*d*m6*sin(d*m6)*besselk(1, a*m6))/((- d^2*m6^2 + 36*pi^2)*(sin(2*m6)/(4*m6) + 1/2)^(1/2)) - (2*a7*d*m7*sin(d*m7)*besselk(1, a*m7))/((- d^2*m7^2 + 36*pi^2)*(sin(2*m7)/(4*m7) + 1/2)^(1/2)) - (2*a8*d*m8*sin(d*m8)*besselk(1, a*m8))/((- d^2*m8^2 + 36*pi^2)*(sin(2*m8)/(4*m8) + 1/2)^(1/2)) - (2*a9*d*m9*sin(d*m9)*besselk(1, a*m9))/((- d^2*m9^2 + 36*pi^2)*(sin(2*m9)/(4*m9) + 1/2)^(1/2)) + (2*a0*d*m0*sinh(d*m0)*(besselj(1, a*m0) + bessely(1, a*m0)*i))/((d^2*m0^2 + 36*pi^2)*(sinh(2*m0)/(4*m0) + 1/2)^(1/2))
 
 
equation7=...
 ...
alpha7 == (2*a1*d*m1*sin(d*m1)*besselk(1, a*m1))/((- d^2*m1^2 + 49*pi^2)*(sin(2*m1)/(4*m1) + 1/2)^(1/2)) + (2*a2*d*m2*sin(d*m2)*besselk(1, a*m2))/((- d^2*m2^2 + 49*pi^2)*(sin(2*m2)/(4*m2) + 1/2)^(1/2)) + (2*a3*d*m3*sin(d*m3)*besselk(1, a*m3))/((- d^2*m3^2 + 49*pi^2)*(sin(2*m3)/(4*m3) + 1/2)^(1/2)) + (2*a4*d*m4*sin(d*m4)*besselk(1, a*m4))/((- d^2*m4^2 + 49*pi^2)*(sin(2*m4)/(4*m4) + 1/2)^(1/2)) + (2*a5*d*m5*sin(d*m5)*besselk(1, a*m5))/((- d^2*m5^2 + 49*pi^2)*(sin(2*m5)/(4*m5) + 1/2)^(1/2)) + (2*a6*d*m6*sin(d*m6)*besselk(1, a*m6))/((- d^2*m6^2 + 49*pi^2)*(sin(2*m6)/(4*m6) + 1/2)^(1/2)) + (2*a7*d*m7*sin(d*m7)*besselk(1, a*m7))/((- d^2*m7^2 + 49*pi^2)*(sin(2*m7)/(4*m7) + 1/2)^(1/2)) + (2*a8*d*m8*sin(d*m8)*besselk(1, a*m8))/((- d^2*m8^2 + 49*pi^2)*(sin(2*m8)/(4*m8) + 1/2)^(1/2)) + (2*a9*d*m9*sin(d*m9)*besselk(1, a*m9))/((- d^2*m9^2 + 49*pi^2)*(sin(2*m9)/(4*m9) + 1/2)^(1/2)) - (2*a0*d*m0*sinh(d*m0)*(besselj(1, a*m0) + bessely(1, a*m0)*i))/((d^2*m0^2 + 49*pi^2)*(sinh(2*m0)/(4*m0) + 1/2)^(1/2))
 
 
equation8=...
 ...
alpha8 == - (2*a1*d*m1*sin(d*m1)*besselk(1, a*m1))/((- d^2*m1^2 + 64*pi^2)*(sin(2*m1)/(4*m1) + 1/2)^(1/2)) - (2*a2*d*m2*sin(d*m2)*besselk(1, a*m2))/((- d^2*m2^2 + 64*pi^2)*(sin(2*m2)/(4*m2) + 1/2)^(1/2)) - (2*a3*d*m3*sin(d*m3)*besselk(1, a*m3))/((- d^2*m3^2 + 64*pi^2)*(sin(2*m3)/(4*m3) + 1/2)^(1/2)) - (2*a4*d*m4*sin(d*m4)*besselk(1, a*m4))/((- d^2*m4^2 + 64*pi^2)*(sin(2*m4)/(4*m4) + 1/2)^(1/2)) - (2*a5*d*m5*sin(d*m5)*besselk(1, a*m5))/((- d^2*m5^2 + 64*pi^2)*(sin(2*m5)/(4*m5) + 1/2)^(1/2)) - (2*a6*d*m6*sin(d*m6)*besselk(1, a*m6))/((- d^2*m6^2 + 64*pi^2)*(sin(2*m6)/(4*m6) + 1/2)^(1/2)) - (2*a7*d*m7*sin(d*m7)*besselk(1, a*m7))/((- d^2*m7^2 + 64*pi^2)*(sin(2*m7)/(4*m7) + 1/2)^(1/2)) - (2*a8*d*m8*sin(d*m8)*besselk(1, a*m8))/((- d^2*m8^2 + 64*pi^2)*(sin(2*m8)/(4*m8) + 1/2)^(1/2)) - (2*a9*d*m9*sin(d*m9)*besselk(1, a*m9))/((- d^2*m9^2 + 64*pi^2)*(sin(2*m9)/(4*m9) + 1/2)^(1/2)) + (2*a0*d*m0*sinh(d*m0)*(besselj(1, a*m0) + bessely(1, a*m0)*i))/((d^2*m0^2 + 64*pi^2)*(sinh(2*m0)/(4*m0) + 1/2)^(1/2))
 
 
equation9=...
 ...
alpha9 == (2*a1*d*m1*sin(d*m1)*besselk(1, a*m1))/((- d^2*m1^2 + 81*pi^2)*(sin(2*m1)/(4*m1) + 1/2)^(1/2)) + (2*a2*d*m2*sin(d*m2)*besselk(1, a*m2))/((- d^2*m2^2 + 81*pi^2)*(sin(2*m2)/(4*m2) + 1/2)^(1/2)) + (2*a3*d*m3*sin(d*m3)*besselk(1, a*m3))/((- d^2*m3^2 + 81*pi^2)*(sin(2*m3)/(4*m3) + 1/2)^(1/2)) + (2*a4*d*m4*sin(d*m4)*besselk(1, a*m4))/((- d^2*m4^2 + 81*pi^2)*(sin(2*m4)/(4*m4) + 1/2)^(1/2)) + (2*a5*d*m5*sin(d*m5)*besselk(1, a*m5))/((- d^2*m5^2 + 81*pi^2)*(sin(2*m5)/(4*m5) + 1/2)^(1/2)) + (2*a6*d*m6*sin(d*m6)*besselk(1, a*m6))/((- d^2*m6^2 + 81*pi^2)*(sin(2*m6)/(4*m6) + 1/2)^(1/2)) + (2*a7*d*m7*sin(d*m7)*besselk(1, a*m7))/((- d^2*m7^2 + 81*pi^2)*(sin(2*m7)/(4*m7) + 1/2)^(1/2)) + (2*a8*d*m8*sin(d*m8)*besselk(1, a*m8))/((- d^2*m8^2 + 81*pi^2)*(sin(2*m8)/(4*m8) + 1/2)^(1/2)) + (2*a9*d*m9*sin(d*m9)*besselk(1, a*m9))/((- d^2*m9^2 + 81*pi^2)*(sin(2*m9)/(4*m9) + 1/2)^(1/2)) - (2*a0*d*m0*sinh(d*m0)*(besselj(1, a*m0) + bessely(1, a*m0)*i))/((d^2*m0^2 + 81*pi^2)*(sinh(2*m0)/(4*m0) + 1/2)^(1/2))
 
 
equation10=...
 ...
a0 == (sinh(sinh(d*m0) - m0)/(m0*(sinh(2*m0)/(4*m0) + 1/2)^(1/2)) - (alpha0*sinh(d*m0))/(2*d*m0*(sinh(2*m0)/(4*m0) + 1/2)^(1/2)) + (pi*alpha1*d*m0*sinh(d*m0)*(besseli(2, m1) + besseli(1, m1)/m1))/((sinh(2*m0)/(4*m0) + 1/2)^(1/2)*besseli(1, m1)*(d^2*m0^2 + pi^2)) - (2*pi*alpha2*d*m0*sinh(d*m0)*(besseli(2, m2) + besseli(1, m2)/m2))/((d^2*m0^2 + 4*pi^2)*(sinh(2*m0)/(4*m0) + 1/2)^(1/2)*besseli(1, m2)) + (3*pi*alpha3*d*m0*sinh(d*m0)*(besseli(2, m3) + besseli(1, m3)/m3))/((d^2*m0^2 + 9*pi^2)*(sinh(2*m0)/(4*m0) + 1/2)^(1/2)*besseli(1, m3)) - (4*pi*alpha4*d*m0*sinh(d*m0)*(besseli(2, m4) + besseli(1, m4)/m4))/((d^2*m0^2 + 16*pi^2)*(sinh(2*m0)/(4*m0) + 1/2)^(1/2)*besseli(1, m4)) + (5*pi*alpha5*d*m0*sinh(d*m0)*(besseli(2, m5) + besseli(1, m5)/m5))/((d^2*m0^2 + 25*pi^2)*(sinh(2*m0)/(4*m0) + 1/2)^(1/2)*besseli(1, m5)) - (6*pi*alpha6*d*m0*sinh(d*m0)*(besseli(2, m6) + besseli(1, m6)/m6))/((d^2*m0^2 + 36*pi^2)*(sinh(2*m0)/(4*m0) + 1/2)^(1/2)*besseli(1, m6)) + (7*pi*alpha7*d*m0*sinh(d*m0)*(besseli(2, m7) + besseli(1, m7)/m7))/((d^2*m0^2 + 49*pi^2)*(sinh(2*m0)/(4*m0) + 1/2)^(1/2)*besseli(1, m7)) - (8*pi*alpha8*d*m0*sinh(d*m0)*(besseli(2, m8) + besseli(1, m8)/m8))/((d^2*m0^2 + 64*pi^2)*(sinh(2*m0)/(4*m0) + 1/2)^(1/2)*besseli(1, m8)) + (9*pi*alpha9*d*m0*sinh(d*m0)*(besseli(2, m9) + besseli(1, m9)/m9))/((d^2*m0^2 + 81*pi^2)*(sinh(2*m0)/(4*m0) + 1/2)^(1/2)*besseli(1, m9)))/(m0*(besselj(2, a*m0) + bessely(2, a*m0)*i - besselj(1, a*m0)/(a*m0) - (bessely(1, a*m0)*i)/(a*m0)))
 
 
equation11=...
 ...
a1 == -((alpha0*sin(d*m1))/(2*d*m1*(sin(2*m1)/(4*m1) + 1/2)^(1/2)) - sin(sin(d*m1) - m1)/(m1*(sin(2*m1)/(4*m1) + 1/2)^(1/2)) + (pi*alpha1*d*m1*sin(d*m1)*(besseli(2, m1) + besseli(1, m1)/m1))/((sin(2*m1)/(4*m1) + 1/2)^(1/2)*besseli(1, m1)*(- d^2*m1^2 + pi^2)) - (2*pi*alpha2*d*m1*sin(d*m1)*(besseli(2, m2) + besseli(1, m2)/m2))/((- d^2*m1^2 + 4*pi^2)*(sin(2*m1)/(4*m1) + 1/2)^(1/2)*besseli(1, m2)) + (3*pi*alpha3*d*m1*sin(d*m1)*(besseli(2, m3) + besseli(1, m3)/m3))/((- d^2*m1^2 + 9*pi^2)*(sin(2*m1)/(4*m1) + 1/2)^(1/2)*besseli(1, m3)) - (4*pi*alpha4*d*m1*sin(d*m1)*(besseli(2, m4) + besseli(1, m4)/m4))/((- d^2*m1^2 + 16*pi^2)*(sin(2*m1)/(4*m1) + 1/2)^(1/2)*besseli(1, m4)) + (5*pi*alpha5*d*m1*sin(d*m1)*(besseli(2, m5) + besseli(1, m5)/m5))/((- d^2*m1^2 + 25*pi^2)*(sin(2*m1)/(4*m1) + 1/2)^(1/2)*besseli(1, m5)) - (6*pi*alpha6*d*m1*sin(d*m1)*(besseli(2, m6) + besseli(1, m6)/m6))/((- d^2*m1^2 + 36*pi^2)*(sin(2*m1)/(4*m1) + 1/2)^(1/2)*besseli(1, m6)) + (7*pi*alpha7*d*m1*sin(d*m1)*(besseli(2, m7) + besseli(1, m7)/m7))/((- d^2*m1^2 + 49*pi^2)*(sin(2*m1)/(4*m1) + 1/2)^(1/2)*besseli(1, m7)) - (8*pi*alpha8*d*m1*sin(d*m1)*(besseli(2, m8) + besseli(1, m8)/m8))/((- d^2*m1^2 + 64*pi^2)*(sin(2*m1)/(4*m1) + 1/2)^(1/2)*besseli(1, m8)) + (9*pi*alpha9*d*m1*sin(d*m1)*(besseli(2, m9) + besseli(1, m9)/m9))/((- d^2*m1^2 + 81*pi^2)*(sin(2*m1)/(4*m1) + 1/2)^(1/2)*besseli(1, m9)))/(m1*(besselk(2, a*m1) - besselk(1, a*m1)/(a*m1)))
 
 
equation12=...
 ...
a2 == -((alpha0*sin(d*m2))/(2*d*m2*(sin(2*m2)/(4*m2) + 1/2)^(1/2)) - sin(sin(d*m2) - m2)/(m2*(sin(2*m2)/(4*m2) + 1/2)^(1/2)) + (pi*alpha1*d*m2*sin(d*m2)*(besseli(2, m1) + besseli(1, m1)/m1))/((sin(2*m2)/(4*m2) + 1/2)^(1/2)*besseli(1, m1)*(- d^2*m2^2 + pi^2)) - (2*pi*alpha2*d*m2*sin(d*m2)*(besseli(2, m2) + besseli(1, m2)/m2))/((- d^2*m2^2 + 4*pi^2)*(sin(2*m2)/(4*m2) + 1/2)^(1/2)*besseli(1, m2)) + (3*pi*alpha3*d*m2*sin(d*m2)*(besseli(2, m3) + besseli(1, m3)/m3))/((- d^2*m2^2 + 9*pi^2)*(sin(2*m2)/(4*m2) + 1/2)^(1/2)*besseli(1, m3)) - (4*pi*alpha4*d*m2*sin(d*m2)*(besseli(2, m4) + besseli(1, m4)/m4))/((- d^2*m2^2 + 16*pi^2)*(sin(2*m2)/(4*m2) + 1/2)^(1/2)*besseli(1, m4)) + (5*pi*alpha5*d*m2*sin(d*m2)*(besseli(2, m5) + besseli(1, m5)/m5))/((- d^2*m2^2 + 25*pi^2)*(sin(2*m2)/(4*m2) + 1/2)^(1/2)*besseli(1, m5)) - (6*pi*alpha6*d*m2*sin(d*m2)*(besseli(2, m6) + besseli(1, m6)/m6))/((- d^2*m2^2 + 36*pi^2)*(sin(2*m2)/(4*m2) + 1/2)^(1/2)*besseli(1, m6)) + (7*pi*alpha7*d*m2*sin(d*m2)*(besseli(2, m7) + besseli(1, m7)/m7))/((- d^2*m2^2 + 49*pi^2)*(sin(2*m2)/(4*m2) + 1/2)^(1/2)*besseli(1, m7)) - (8*pi*alpha8*d*m2*sin(d*m2)*(besseli(2, m8) + besseli(1, m8)/m8))/((- d^2*m2^2 + 64*pi^2)*(sin(2*m2)/(4*m2) + 1/2)^(1/2)*besseli(1, m8)) + (9*pi*alpha9*d*m2*sin(d*m2)*(besseli(2, m9) + besseli(1, m9)/m9))/((- d^2*m2^2 + 81*pi^2)*(sin(2*m2)/(4*m2) + 1/2)^(1/2)*besseli(1, m9)))/(m2*(besselk(2, a*m2) - besselk(1, a*m2)/(a*m2)))
 
 
equation13=...
 ...
a3 == -((alpha0*sin(d*m3))/(2*d*m3*(sin(2*m3)/(4*m3) + 1/2)^(1/2)) - sin(sin(d*m3) - m3)/(m3*(sin(2*m3)/(4*m3) + 1/2)^(1/2)) + (pi*alpha1*d*m3*sin(d*m3)*(besseli(2, m1) + besseli(1, m1)/m1))/((sin(2*m3)/(4*m3) + 1/2)^(1/2)*besseli(1, m1)*(- d^2*m3^2 + pi^2)) - (2*pi*alpha2*d*m3*sin(d*m3)*(besseli(2, m2) + besseli(1, m2)/m2))/((- d^2*m3^2 + 4*pi^2)*(sin(2*m3)/(4*m3) + 1/2)^(1/2)*besseli(1, m2)) + (3*pi*alpha3*d*m3*sin(d*m3)*(besseli(2, m3) + besseli(1, m3)/m3))/((- d^2*m3^2 + 9*pi^2)*(sin(2*m3)/(4*m3) + 1/2)^(1/2)*besseli(1, m3)) - (4*pi*alpha4*d*m3*sin(d*m3)*(besseli(2, m4) + besseli(1, m4)/m4))/((- d^2*m3^2 + 16*pi^2)*(sin(2*m3)/(4*m3) + 1/2)^(1/2)*besseli(1, m4)) + (5*pi*alpha5*d*m3*sin(d*m3)*(besseli(2, m5) + besseli(1, m5)/m5))/((- d^2*m3^2 + 25*pi^2)*(sin(2*m3)/(4*m3) + 1/2)^(1/2)*besseli(1, m5)) - (6*pi*alpha6*d*m3*sin(d*m3)*(besseli(2, m6) + besseli(1, m6)/m6))/((- d^2*m3^2 + 36*pi^2)*(sin(2*m3)/(4*m3) + 1/2)^(1/2)*besseli(1, m6)) + (7*pi*alpha7*d*m3*sin(d*m3)*(besseli(2, m7) + besseli(1, m7)/m7))/((- d^2*m3^2 + 49*pi^2)*(sin(2*m3)/(4*m3) + 1/2)^(1/2)*besseli(1, m7)) - (8*pi*alpha8*d*m3*sin(d*m3)*(besseli(2, m8) + besseli(1, m8)/m8))/((- d^2*m3^2 + 64*pi^2)*(sin(2*m3)/(4*m3) + 1/2)^(1/2)*besseli(1, m8)) + (9*pi*alpha9*d*m3*sin(d*m3)*(besseli(2, m9) + besseli(1, m9)/m9))/((- d^2*m3^2 + 81*pi^2)*(sin(2*m3)/(4*m3) + 1/2)^(1/2)*besseli(1, m9)))/(m3*(besselk(2, a*m3) - besselk(1, a*m3)/(a*m3)))
 
 
equation14=...
 ...
a4 == -((alpha0*sin(d*m4))/(2*d*m4*(sin(2*m4)/(4*m4) + 1/2)^(1/2)) - sin(sin(d*m4) - m4)/(m4*(sin(2*m4)/(4*m4) + 1/2)^(1/2)) + (pi*alpha1*d*m4*sin(d*m4)*(besseli(2, m1) + besseli(1, m1)/m1))/((sin(2*m4)/(4*m4) + 1/2)^(1/2)*besseli(1, m1)*(- d^2*m4^2 + pi^2)) - (2*pi*alpha2*d*m4*sin(d*m4)*(besseli(2, m2) + besseli(1, m2)/m2))/((- d^2*m4^2 + 4*pi^2)*(sin(2*m4)/(4*m4) + 1/2)^(1/2)*besseli(1, m2)) + (3*pi*alpha3*d*m4*sin(d*m4)*(besseli(2, m3) + besseli(1, m3)/m3))/((- d^2*m4^2 + 9*pi^2)*(sin(2*m4)/(4*m4) + 1/2)^(1/2)*besseli(1, m3)) - (4*pi*alpha4*d*m4*sin(d*m4)*(besseli(2, m4) + besseli(1, m4)/m4))/((- d^2*m4^2 + 16*pi^2)*(sin(2*m4)/(4*m4) + 1/2)^(1/2)*besseli(1, m4)) + (5*pi*alpha5*d*m4*sin(d*m4)*(besseli(2, m5) + besseli(1, m5)/m5))/((- d^2*m4^2 + 25*pi^2)*(sin(2*m4)/(4*m4) + 1/2)^(1/2)*besseli(1, m5)) - (6*pi*alpha6*d*m4*sin(d*m4)*(besseli(2, m6) + besseli(1, m6)/m6))/((- d^2*m4^2 + 36*pi^2)*(sin(2*m4)/(4*m4) + 1/2)^(1/2)*besseli(1, m6)) + (7*pi*alpha7*d*m4*sin(d*m4)*(besseli(2, m7) + besseli(1, m7)/m7))/((- d^2*m4^2 + 49*pi^2)*(sin(2*m4)/(4*m4) + 1/2)^(1/2)*besseli(1, m7)) - (8*pi*alpha8*d*m4*sin(d*m4)*(besseli(2, m8) + besseli(1, m8)/m8))/((- d^2*m4^2 + 64*pi^2)*(sin(2*m4)/(4*m4) + 1/2)^(1/2)*besseli(1, m8)) + (9*pi*alpha9*d*m4*sin(d*m4)*(besseli(2, m9) + besseli(1, m9)/m9))/((- d^2*m4^2 + 81*pi^2)*(sin(2*m4)/(4*m4) + 1/2)^(1/2)*besseli(1, m9)))/(m4*(besselk(2, a*m4) - besselk(1, a*m4)/(a*m4)))
 
 
equation15=...
 ...
a5 == -((alpha0*sin(d*m5))/(2*d*m5*(sin(2*m5)/(4*m5) + 1/2)^(1/2)) - sin(sin(d*m5) - m5)/(m5*(sin(2*m5)/(4*m5) + 1/2)^(1/2)) + (pi*alpha1*d*m5*sin(d*m5)*(besseli(2, m1) + besseli(1, m1)/m1))/((sin(2*m5)/(4*m5) + 1/2)^(1/2)*besseli(1, m1)*(- d^2*m5^2 + pi^2)) - (2*pi*alpha2*d*m5*sin(d*m5)*(besseli(2, m2) + besseli(1, m2)/m2))/((- d^2*m5^2 + 4*pi^2)*(sin(2*m5)/(4*m5) + 1/2)^(1/2)*besseli(1, m2)) + (3*pi*alpha3*d*m5*sin(d*m5)*(besseli(2, m3) + besseli(1, m3)/m3))/((- d^2*m5^2 + 9*pi^2)*(sin(2*m5)/(4*m5) + 1/2)^(1/2)*besseli(1, m3)) - (4*pi*alpha4*d*m5*sin(d*m5)*(besseli(2, m4) + besseli(1, m4)/m4))/((- d^2*m5^2 + 16*pi^2)*(sin(2*m5)/(4*m5) + 1/2)^(1/2)*besseli(1, m4)) + (5*pi*alpha5*d*m5*sin(d*m5)*(besseli(2, m5) + besseli(1, m5)/m5))/((- d^2*m5^2 + 25*pi^2)*(sin(2*m5)/(4*m5) + 1/2)^(1/2)*besseli(1, m5)) - (6*pi*alpha6*d*m5*sin(d*m5)*(besseli(2, m6) + besseli(1, m6)/m6))/((- d^2*m5^2 + 36*pi^2)*(sin(2*m5)/(4*m5) + 1/2)^(1/2)*besseli(1, m6)) + (7*pi*alpha7*d*m5*sin(d*m5)*(besseli(2, m7) + besseli(1, m7)/m7))/((- d^2*m5^2 + 49*pi^2)*(sin(2*m5)/(4*m5) + 1/2)^(1/2)*besseli(1, m7)) - (8*pi*alpha8*d*m5*sin(d*m5)*(besseli(2, m8) + besseli(1, m8)/m8))/((- d^2*m5^2 + 64*pi^2)*(sin(2*m5)/(4*m5) + 1/2)^(1/2)*besseli(1, m8)) + (9*pi*alpha9*d*m5*sin(d*m5)*(besseli(2, m9) + besseli(1, m9)/m9))/((- d^2*m5^2 + 81*pi^2)*(sin(2*m5)/(4*m5) + 1/2)^(1/2)*besseli(1, m9)))/(m5*(besselk(2, a*m5) - besselk(1, a*m5)/(a*m5)))
 
 
equation16=...
 ...
a6 == -((alpha0*sin(d*m6))/(2*d*m6*(sin(2*m6)/(4*m6) + 1/2)^(1/2)) - sin(sin(d*m6) - m6)/(m6*(sin(2*m6)/(4*m6) + 1/2)^(1/2)) + (pi*alpha1*d*m6*sin(d*m6)*(besseli(2, m1) + besseli(1, m1)/m1))/((sin(2*m6)/(4*m6) + 1/2)^(1/2)*besseli(1, m1)*(- d^2*m6^2 + pi^2)) - (2*pi*alpha2*d*m6*sin(d*m6)*(besseli(2, m2) + besseli(1, m2)/m2))/((- d^2*m6^2 + 4*pi^2)*(sin(2*m6)/(4*m6) + 1/2)^(1/2)*besseli(1, m2)) + (3*pi*alpha3*d*m6*sin(d*m6)*(besseli(2, m3) + besseli(1, m3)/m3))/((- d^2*m6^2 + 9*pi^2)*(sin(2*m6)/(4*m6) + 1/2)^(1/2)*besseli(1, m3)) - (4*pi*alpha4*d*m6*sin(d*m6)*(besseli(2, m4) + besseli(1, m4)/m4))/((- d^2*m6^2 + 16*pi^2)*(sin(2*m6)/(4*m6) + 1/2)^(1/2)*besseli(1, m4)) + (5*pi*alpha5*d*m6*sin(d*m6)*(besseli(2, m5) + besseli(1, m5)/m5))/((- d^2*m6^2 + 25*pi^2)*(sin(2*m6)/(4*m6) + 1/2)^(1/2)*besseli(1, m5)) - (6*pi*alpha6*d*m6*sin(d*m6)*(besseli(2, m6) + besseli(1, m6)/m6))/((- d^2*m6^2 + 36*pi^2)*(sin(2*m6)/(4*m6) + 1/2)^(1/2)*besseli(1, m6)) + (7*pi*alpha7*d*m6*sin(d*m6)*(besseli(2, m7) + besseli(1, m7)/m7))/((- d^2*m6^2 + 49*pi^2)*(sin(2*m6)/(4*m6) + 1/2)^(1/2)*besseli(1, m7)) - (8*pi*alpha8*d*m6*sin(d*m6)*(besseli(2, m8) + besseli(1, m8)/m8))/((- d^2*m6^2 + 64*pi^2)*(sin(2*m6)/(4*m6) + 1/2)^(1/2)*besseli(1, m8)) + (9*pi*alpha9*d*m6*sin(d*m6)*(besseli(2, m9) + besseli(1, m9)/m9))/((- d^2*m6^2 + 81*pi^2)*(sin(2*m6)/(4*m6) + 1/2)^(1/2)*besseli(1, m9)))/(m6*(besselk(2, a*m6) - besselk(1, a*m6)/(a*m6)))
 
 
equation17=...
 ...
a7 == -((alpha0*sin(d*m7))/(2*d*m7*(sin(2*m7)/(4*m7) + 1/2)^(1/2)) - sin(sin(d*m7) - m7)/(m7*(sin(2*m7)/(4*m7) + 1/2)^(1/2)) + (pi*alpha1*d*m7*sin(d*m7)*(besseli(2, m1) + besseli(1, m1)/m1))/((sin(2*m7)/(4*m7) + 1/2)^(1/2)*besseli(1, m1)*(- d^2*m7^2 + pi^2)) - (2*pi*alpha2*d*m7*sin(d*m7)*(besseli(2, m2) + besseli(1, m2)/m2))/((- d^2*m7^2 + 4*pi^2)*(sin(2*m7)/(4*m7) + 1/2)^(1/2)*besseli(1, m2)) + (3*pi*alpha3*d*m7*sin(d*m7)*(besseli(2, m3) + besseli(1, m3)/m3))/((- d^2*m7^2 + 9*pi^2)*(sin(2*m7)/(4*m7) + 1/2)^(1/2)*besseli(1, m3)) - (4*pi*alpha4*d*m7*sin(d*m7)*(besseli(2, m4) + besseli(1, m4)/m4))/((- d^2*m7^2 + 16*pi^2)*(sin(2*m7)/(4*m7) + 1/2)^(1/2)*besseli(1, m4)) + (5*pi*alpha5*d*m7*sin(d*m7)*(besseli(2, m5) + besseli(1, m5)/m5))/((- d^2*m7^2 + 25*pi^2)*(sin(2*m7)/(4*m7) + 1/2)^(1/2)*besseli(1, m5)) - (6*pi*alpha6*d*m7*sin(d*m7)*(besseli(2, m6) + besseli(1, m6)/m6))/((- d^2*m7^2 + 36*pi^2)*(sin(2*m7)/(4*m7) + 1/2)^(1/2)*besseli(1, m6)) + (7*pi*alpha7*d*m7*sin(d*m7)*(besseli(2, m7) + besseli(1, m7)/m7))/((- d^2*m7^2 + 49*pi^2)*(sin(2*m7)/(4*m7) + 1/2)^(1/2)*besseli(1, m7)) - (8*pi*alpha8*d*m7*sin(d*m7)*(besseli(2, m8) + besseli(1, m8)/m8))/((- d^2*m7^2 + 64*pi^2)*(sin(2*m7)/(4*m7) + 1/2)^(1/2)*besseli(1, m8)) + (9*pi*alpha9*d*m7*sin(d*m7)*(besseli(2, m9) + besseli(1, m9)/m9))/((- d^2*m7^2 + 81*pi^2)*(sin(2*m7)/(4*m7) + 1/2)^(1/2)*besseli(1, m9)))/(m7*(besselk(2, a*m7) - besselk(1, a*m7)/(a*m7)))
 
 
equation18=...
 ...
a8 == -((alpha0*sin(d*m8))/(2*d*m8*(sin(2*m8)/(4*m8) + 1/2)^(1/2)) - sin(sin(d*m8) - m8)/(m8*(sin(2*m8)/(4*m8) + 1/2)^(1/2)) + (pi*alpha1*d*m8*sin(d*m8)*(besseli(2, m1) + besseli(1, m1)/m1))/((sin(2*m8)/(4*m8) + 1/2)^(1/2)*besseli(1, m1)*(- d^2*m8^2 + pi^2)) - (2*pi*alpha2*d*m8*sin(d*m8)*(besseli(2, m2) + besseli(1, m2)/m2))/((- d^2*m8^2 + 4*pi^2)*(sin(2*m8)/(4*m8) + 1/2)^(1/2)*besseli(1, m2)) + (3*pi*alpha3*d*m8*sin(d*m8)*(besseli(2, m3) + besseli(1, m3)/m3))/((- d^2*m8^2 + 9*pi^2)*(sin(2*m8)/(4*m8) + 1/2)^(1/2)*besseli(1, m3)) - (4*pi*alpha4*d*m8*sin(d*m8)*(besseli(2, m4) + besseli(1, m4)/m4))/((- d^2*m8^2 + 16*pi^2)*(sin(2*m8)/(4*m8) + 1/2)^(1/2)*besseli(1, m4)) + (5*pi*alpha5*d*m8*sin(d*m8)*(besseli(2, m5) + besseli(1, m5)/m5))/((- d^2*m8^2 + 25*pi^2)*(sin(2*m8)/(4*m8) + 1/2)^(1/2)*besseli(1, m5)) - (6*pi*alpha6*d*m8*sin(d*m8)*(besseli(2, m6) + besseli(1, m6)/m6))/((- d^2*m8^2 + 36*pi^2)*(sin(2*m8)/(4*m8) + 1/2)^(1/2)*besseli(1, m6)) + (7*pi*alpha7*d*m8*sin(d*m8)*(besseli(2, m7) + besseli(1, m7)/m7))/((- d^2*m8^2 + 49*pi^2)*(sin(2*m8)/(4*m8) + 1/2)^(1/2)*besseli(1, m7)) - (8*pi*alpha8*d*m8*sin(d*m8)*(besseli(2, m8) + besseli(1, m8)/m8))/((- d^2*m8^2 + 64*pi^2)*(sin(2*m8)/(4*m8) + 1/2)^(1/2)*besseli(1, m8)) + (9*pi*alpha9*d*m8*sin(d*m8)*(besseli(2, m9) + besseli(1, m9)/m9))/((- d^2*m8^2 + 81*pi^2)*(sin(2*m8)/(4*m8) + 1/2)^(1/2)*besseli(1, m9)))/(m8*(besselk(2, a*m8) - besselk(1, a*m8)/(a*m8)))
 
 
equation19=...
 ...
a9 == -((alpha0*sin(d*m9))/(2*d*m9*(sin(2*m9)/(4*m9) + 1/2)^(1/2)) - sin(sin(d*m9) - m9)/(m9*(sin(2*m9)/(4*m9) + 1/2)^(1/2)) + (pi*alpha1*d*m9*sin(d*m9)*(besseli(2, m1) + besseli(1, m1)/m1))/((sin(2*m9)/(4*m9) + 1/2)^(1/2)*besseli(1, m1)*(- d^2*m9^2 + pi^2)) - (2*pi*alpha2*d*m9*sin(d*m9)*(besseli(2, m2) + besseli(1, m2)/m2))/((- d^2*m9^2 + 4*pi^2)*(sin(2*m9)/(4*m9) + 1/2)^(1/2)*besseli(1, m2)) + (3*pi*alpha3*d*m9*sin(d*m9)*(besseli(2, m3) + besseli(1, m3)/m3))/((- d^2*m9^2 + 9*pi^2)*(sin(2*m9)/(4*m9) + 1/2)^(1/2)*besseli(1, m3)) - (4*pi*alpha4*d*m9*sin(d*m9)*(besseli(2, m4) + besseli(1, m4)/m4))/((- d^2*m9^2 + 16*pi^2)*(sin(2*m9)/(4*m9) + 1/2)^(1/2)*besseli(1, m4)) + (5*pi*alpha5*d*m9*sin(d*m9)*(besseli(2, m5) + besseli(1, m5)/m5))/((- d^2*m9^2 + 25*pi^2)*(sin(2*m9)/(4*m9) + 1/2)^(1/2)*besseli(1, m5)) - (6*pi*alpha6*d*m9*sin(d*m9)*(besseli(2, m6) + besseli(1, m6)/m6))/((- d^2*m9^2 + 36*pi^2)*(sin(2*m9)/(4*m9) + 1/2)^(1/2)*besseli(1, m6)) + (7*pi*alpha7*d*m9*sin(d*m9)*(besseli(2, m7) + besseli(1, m7)/m7))/((- d^2*m9^2 + 49*pi^2)*(sin(2*m9)/(4*m9) + 1/2)^(1/2)*besseli(1, m7)) - (8*pi*alpha8*d*m9*sin(d*m9)*(besseli(2, m8) + besseli(1, m8)/m8))/((- d^2*m9^2 + 64*pi^2)*(sin(2*m9)/(4*m9) + 1/2)^(1/2)*besseli(1, m8)) + (9*pi*alpha9*d*m9*sin(d*m9)*(besseli(2, m9) + besseli(1, m9)/m9))/((- d^2*m9^2 + 81*pi^2)*(sin(2*m9)/(4*m9) + 1/2)^(1/2)*besseli(1, m9)))/(m9*(besselk(2, a*m9) - besselk(1, a*m9)/(a*m9)))
 
 
[A, b] = equationsToMatrix([equation0, equation1, equation2,equation3,equation4,equation5,equation6, equation7, equation8,equation9,equation10,equation11,equation12, equation13, equation14,equation15,equation16,equation17,equation18, equation19], [a0, a1, a2 ,a3, a4 ,a5,a6,a7,a8,a9, alpha0 ,alpha1 ,alpha2 ,alpha3 ,alpha4, alpha5,alpha6 ,alpha7 ,alpha8, alpha9]);
 
x = A\b
 

[dcydhb dcydhb]

and why when i use 'vpasolve',it runs out the root that

Unknown quantity contains other unknown quantity

for example

[a0, a1,alpha0 ,alpha1]=vpasolve([equation0, equation1, equation2,equation3], [a0, a1,alpha0 ,alpha1])

[Jan]

I'm not sure how smart Matlab handles the power operation. But in general ^(1/2) is much more expensive than sqrt().

[dcydhb dcydhb]

I want to know how to speed up structurally or where it affects speed the most.

[Jan]

Then use the profile command, which tells you the details exactly.