omega=1.795200;
m0=7.234737;
m1=1.816835;
m2=5.349043;
m3=8.732933;
m4=12.024736;
m5=15.265388;
m6=18.476340;
m7=21.668910;
m8=24.849430;
m9=28.021668;
m10=31.187986;
syms a0 a1 a2 a3 a4 a5 a6 a7 a8 a9 alpha0 alpha1 alpha2 alpha3 alpha4 alpha5 alpha6 alpha7 alpha8 alpha9;
h=22;
hgang=22;
agang=2.2497;
dgang=21.1940;
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*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)*(pi^2 - d^2*m1^2)) + (2*a2*d*m2*sin(d*m2)*besselk(1, a*m2))/((sin(2*m2)/(4*m2) + 1/2)^(1/2)*(pi^2 - d^2*m2^2)) + (2*a3*d*m3*sin(d*m3)*besselk(1, a*m3))/((sin(2*m3)/(4*m3) + 1/2)^(1/2)*(pi^2 - d^2*m3^2)) + (2*a4*d*m4*sin(d*m4)*besselk(1, a*m4))/((sin(2*m4)/(4*m4) + 1/2)^(1/2)*(pi^2 - d^2*m4^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)*(pi^2 + d^2*m0^2));
equation2=alpha2 == (2*a0*d*m0*sinh(d*m0)*(besselj(1, a*m0) + bessely(1, a*m0)*i))/((4*pi^2 + d^2*m0^2)*(sinh(2*m0)/(4*m0) + 1/2)^(1/2)) - (2*a2*d*m2*sin(d*m2)*besselk(1, a*m2))/((4*pi^2 - d^2*m2^2)*(sin(2*m2)/(4*m2) + 1/2)^(1/2)) - (2*a3*d*m3*sin(d*m3)*besselk(1, a*m3))/((4*pi^2 - d^2*m3^2)*(sin(2*m3)/(4*m3) + 1/2)^(1/2)) - (2*a4*d*m4*sin(d*m4)*besselk(1, a*m4))/((4*pi^2 - d^2*m4^2)*(sin(2*m4)/(4*m4) + 1/2)^(1/2)) - (2*a1*d*m1*sin(d*m1)*besselk(1, a*m1))/((4*pi^2 - d^2*m1^2)*(sin(2*m1)/(4*m1) + 1/2)^(1/2));
equation3=alpha3 == (2*a1*d*m1*sin(d*m1)*besselk(1, a*m1))/((9*pi^2 - d^2*m1^2)*(sin(2*m1)/(4*m1) + 1/2)^(1/2)) + (2*a2*d*m2*sin(d*m2)*besselk(1, a*m2))/((9*pi^2 - d^2*m2^2)*(sin(2*m2)/(4*m2) + 1/2)^(1/2)) + (2*a3*d*m3*sin(d*m3)*besselk(1, a*m3))/((9*pi^2 - d^2*m3^2)*(sin(2*m3)/(4*m3) + 1/2)^(1/2)) + (2*a4*d*m4*sin(d*m4)*besselk(1, a*m4))/((9*pi^2 - d^2*m4^2)*(sin(2*m4)/(4*m4) + 1/2)^(1/2)) - (2*a0*d*m0*sinh(d*m0)*(besselj(1, a*m0) + bessely(1, a*m0)*i))/((9*pi^2 + d^2*m0^2)*(sinh(2*m0)/(4*m0) + 1/2)^(1/2));
equation4=alpha4 == (2*a0*d*m0*sinh(d*m0)*(besselj(1, a*m0) + bessely(1, a*m0)*i))/((16*pi^2 + d^2*m0^2)*(sinh(2*m0)/(4*m0) + 1/2)^(1/2)) - (2*a2*d*m2*sin(d*m2)*besselk(1, a*m2))/((16*pi^2 - d^2*m2^2)*(sin(2*m2)/(4*m2) + 1/2)^(1/2)) - (2*a3*d*m3*sin(d*m3)*besselk(1, a*m3))/((16*pi^2 - d^2*m3^2)*(sin(2*m3)/(4*m3) + 1/2)^(1/2)) - (2*a4*d*m4*sin(d*m4)*besselk(1, a*m4))/((16*pi^2 - d^2*m4^2)*(sin(2*m4)/(4*m4) + 1/2)^(1/2)) - (2*a1*d*m1*sin(d*m1)*besselk(1, a*m1))/((16*pi^2 - d^2*m1^2)*(sin(2*m1)/(4*m1) + 1/2)^(1/2));
equation5=a0 == (sinh(sinh(d*m0) - m0)/(m0*(sinh(2*m0)/(4*m0) + 1/2)^(1/2)) - (alpha0*sinh(d*m0))/(2*a*m0*(sinh(2*m0)/(4*m0) + 1/2)^(1/2)) + (pi*alpha1*d*m0*sinh(d*m0)*(besseli(2, a*m1) + besseli(1, a*m1)/(a*m1)))/(besseli(1, a*m1)*(sinh(2*m0)/(4*m0) + 1/2)^(1/2)*(pi^2 + d^2*m0^2)) - (2*pi*alpha2*d*m0*sinh(d*m0)*(besseli(2, a*m2) + besseli(1, a*m2)/(a*m2)))/((4*pi^2 + d^2*m0^2)*besseli(1, a*m2)*(sinh(2*m0)/(4*m0) + 1/2)^(1/2)) + (3*pi*alpha3*d*m0*sinh(d*m0)*(besseli(2, a*m3) + besseli(1, a*m3)/(a*m3)))/((9*pi^2 + d^2*m0^2)*besseli(1, a*m3)*(sinh(2*m0)/(4*m0) + 1/2)^(1/2)) - (4*pi*alpha4*d*m0*sinh(d*m0)*(besseli(2, a*m4) + besseli(1, a*m4)/(a*m4)))/((16*pi^2 + d^2*m0^2)*besseli(1, a*m4)*(sinh(2*m0)/(4*m0) + 1/2)^(1/2)))/(m0*(besselj(2, a*m0) + bessely(2, a*m0)*i - besselj(1, a*m0)/(a*m0) + (bessely(1, a*m0)*(-i))/(a*m0)));
equation6=a1 == (sin(sin(d*m1) - m1)/(m1*(sin(2*m1)/(4*m1) + 1/2)^(1/2)) - (alpha0*sin(d*m1))/(2*a*m1*(sin(2*m1)/(4*m1) + 1/2)^(1/2)) - (pi*alpha1*d*m1*sin(d*m1)*(besseli(2, a*m1) + besseli(1, a*m1)/(a*m1)))/(besseli(1, a*m1)*(sin(2*m1)/(4*m1) + 1/2)^(1/2)*(pi^2 - d^2*m1^2)) + (2*pi*alpha2*d*m1*sin(d*m1)*(besseli(2, a*m2) + besseli(1, a*m2)/(a*m2)))/((4*pi^2 - d^2*m1^2)*besseli(1, a*m2)*(sin(2*m1)/(4*m1) + 1/2)^(1/2)) - (3*pi*alpha3*d*m1*sin(d*m1)*(besseli(2, a*m3) + besseli(1, a*m3)/(a*m3)))/((9*pi^2 - d^2*m1^2)*besseli(1, a*m3)*(sin(2*m1)/(4*m1) + 1/2)^(1/2)) + (4*pi*alpha4*d*m1*sin(d*m1)*(besseli(2, a*m4) + besseli(1, a*m4)/(a*m4)))/((16*pi^2 - d^2*m1^2)*besseli(1, a*m4)*(sin(2*m1)/(4*m1) + 1/2)^(1/2)))/(m1*(besselk(2, a*m1) - besselk(1, a*m1)/(a*m1)));
equation7=a2 == (sin(sin(d*m2) - m2)/(m2*(sin(2*m2)/(4*m2) + 1/2)^(1/2)) - (alpha0*sin(d*m2))/(2*a*m2*(sin(2*m2)/(4*m2) + 1/2)^(1/2)) - (pi*alpha1*d*m2*sin(d*m2)*(besseli(2, a*m1) + besseli(1, a*m1)/(a*m1)))/(besseli(1, a*m1)*(sin(2*m2)/(4*m2) + 1/2)^(1/2)*(pi^2 - d^2*m2^2)) + (2*pi*alpha2*d*m2*sin(d*m2)*(besseli(2, a*m2) + besseli(1, a*m2)/(a*m2)))/((4*pi^2 - d^2*m2^2)*besseli(1, a*m2)*(sin(2*m2)/(4*m2) + 1/2)^(1/2)) - (3*pi*alpha3*d*m2*sin(d*m2)*(besseli(2, a*m3) + besseli(1, a*m3)/(a*m3)))/((9*pi^2 - d^2*m2^2)*besseli(1, a*m3)*(sin(2*m2)/(4*m2) + 1/2)^(1/2)) + (4*pi*alpha4*d*m2*sin(d*m2)*(besseli(2, a*m4) + besseli(1, a*m4)/(a*m4)))/((16*pi^2 - d^2*m2^2)*besseli(1, a*m4)*(sin(2*m2)/(4*m2) + 1/2)^(1/2)))/(m2*(besselk(2, a*m2) - besselk(1, a*m2)/(a*m2)));
equation8=a3 == (sin(sin(d*m3) - m3)/(m3*(sin(2*m3)/(4*m3) + 1/2)^(1/2)) - (alpha0*sin(d*m3))/(2*a*m3*(sin(2*m3)/(4*m3) + 1/2)^(1/2)) - (pi*alpha1*d*m3*sin(d*m3)*(besseli(2, a*m1) + besseli(1, a*m1)/(a*m1)))/(besseli(1, a*m1)*(sin(2*m3)/(4*m3) + 1/2)^(1/2)*(pi^2 - d^2*m3^2)) + (2*pi*alpha2*d*m3*sin(d*m3)*(besseli(2, a*m2) + besseli(1, a*m2)/(a*m2)))/((4*pi^2 - d^2*m3^2)*besseli(1, a*m2)*(sin(2*m3)/(4*m3) + 1/2)^(1/2)) - (3*pi*alpha3*d*m3*sin(d*m3)*(besseli(2, a*m3) + besseli(1, a*m3)/(a*m3)))/((9*pi^2 - d^2*m3^2)*besseli(1, a*m3)*(sin(2*m3)/(4*m3) + 1/2)^(1/2)) + (4*pi*alpha4*d*m3*sin(d*m3)*(besseli(2, a*m4) + besseli(1, a*m4)/(a*m4)))/((16*pi^2 - d^2*m3^2)*besseli(1, a*m4)*(sin(2*m3)/(4*m3) + 1/2)^(1/2)))/(m3*(besselk(2, a*m3) - besselk(1, a*m3)/(a*m3)));
equation9=a4 == (sin(sin(d*m4) - m4)/(m4*(sin(2*m4)/(4*m4) + 1/2)^(1/2)) - (alpha0*sin(d*m4))/(2*a*m4*(sin(2*m4)/(4*m4) + 1/2)^(1/2)) - (pi*alpha1*d*m4*sin(d*m4)*(besseli(2, a*m1) + besseli(1, a*m1)/(a*m1)))/(besseli(1, a*m1)*(sin(2*m4)/(4*m4) + 1/2)^(1/2)*(pi^2 - d^2*m4^2)) + (2*pi*alpha2*d*m4*sin(d*m4)*(besseli(2, a*m2) + besseli(1, a*m2)/(a*m2)))/((4*pi^2 - d^2*m4^2)*besseli(1, a*m2)*(sin(2*m4)/(4*m4) + 1/2)^(1/2)) - (3*pi*alpha3*d*m4*sin(d*m4)*(besseli(2, a*m3) + besseli(1, a*m3)/(a*m3)))/((9*pi^2 - d^2*m4^2)*besseli(1, a*m3)*(sin(2*m4)/(4*m4) + 1/2)^(1/2)) + (4*pi*alpha4*d*m4*sin(d*m4)*(besseli(2, a*m4) + besseli(1, a*m4)/(a*m4)))/((16*pi^2 - d^2*m4^2)*besseli(1, a*m4)*(sin(2*m4)/(4*m4) + 1/2)^(1/2)))/(m4*(besselk(2, a*m4) - besselk(1, a*m4)/(a*m4)));
[A, b] = equationsToMatrix([equation0, equation1, equation2,equation3,equation4,equation5,equation6, equation7, equation8,equation9], [a0, a1, a2 ,a3, a4 , alpha0 ,alpha1 ,alpha2 ,alpha3 ,alpha4]);
x = A\b;
simplify( A*x-b)
format long
solution=@(n)(vpa(x(n+1),10));
a0= solution(0)
a1= solution(1)
a2 = solution(2)
a3= solution(3)
a4 = solution(4)
alpha0= solution(5)
alpha1 = solution(6)
alpha2 = solution(7)
alpha3 = solution(8)
alpha4 = solution(9)