function An = An2(n)
global a b kp_a kp_b ks_a ks_b mu_a lambda_a mu_b lambda_b omega vs vp rho_b

Cn(:,1) = [kp_b*(-besselh(n+1,1,kp_b*a)+n/(kp_b*a)*besselh(n,1,kp_b*a))
           -n/a*besselh(n,1,kp_b*a)
           kp_b**2*((-lambda_b-2*mu_b/(kp_b*a)**2*(n+(kp_b*a)**2-n**2))*besselh(n,1,kp_b*a)+2*mu_b/(kp_b*a)*besselh(n+1,1,kp_b*a))
           -2*mu_b*n/a*(-kp_b*besselh(n+1,1,kp_b*a)+(n-1)/a*besselh(n,1,kp_b*a))
           2*mu_b/b**2*((n**2-n-0.5*(ks_b*b)**2)*besselh(n,1,kp_b*b)+kp_b*b*besselh(n+1,1,kp_b*b))-1i*omega*rho_b*vp*kp_b*(-besselh(n+1,1,kp_b*b)+n/(kp_b*b)*besselh(n,1,kp_b*b))
           -2*mu_b*n/(b**2)*((n-1)*besselh(n,1,kp_b*b)-kp_b*b*besselh(n+1,1,kp_b*b))+1i*omega*rho_b*vs*n/b*besselh(n,1,kp_b*b)];
           
Cn(:,2) = [n/a*besselh(n,1,ks_b*a)
           -ks_b*(-besselh(n+1,1,ks_b*a)+n/(ks_b*a)*besselh(n,1,ks_b*a))
           2*mu_b*n/a*(-1/a*(1-n)*besselh(n,1,ks_b*a)-ks_b*besselh(n+1,1,ks_b*a))
           mu_b*2/a*(1/a*(-n**2+n+0.5*(ks_b*a)**2)*besselh(n,1,ks_b*a)-ks_b*besselh(n+1,1,ks_b*a))
           2*mu_b*n/b**2*((n-1)*besselh(n,1,ks_b*b)-ks_b*b*besselh(n+1,1,ks_b*b))-1i*omega*rho_b*vp*n/b*besselh(n,1,ks_b*b)
           -2*mu_b/(b**2)*((n**2-n-0.5*(ks_b*b)**2)*besselh(n,1,ks_b*b)+ks_b*b*besselh(n+1,1,ks_b*b))+1i*omega*rho_b*vs*ks_b*(-besselh(n+1,1,ks_b*b)+n/(ks_b*b)*besselh(n,1,ks_b*b))];  
           
           
Cn(:,3)= [kp_b*(besselh(n-1,2,kp_b*a)-n/(kp_b*a)*besselh(n,2,kp_b*a))
          -n/a*besselh(n,2,kp_b*a)
          kp_b**2*((-lambda_b+2*mu_b/(kp_b*a)**2*(n-(kp_b*a)**2+n**2))*besselh(n,2,kp_b*a)-2*mu_b/(kp_b*a)*besselh(n-1,2,kp_b*a))
          -2*mu_b*n/a*(kp_b*besselh(n-1,2,kp_b*a)-(n+1)/a*besselh(n,2,kp_b*a))
          -mu_b*((ks_b**2-2*kp_b**2)*besselh(n,2,kp_b*b)-2*kp_b**2*(besselh(n-2,2,kp_b*b)-(n-1)/(kp_b*b)*besselh(n-1,2,kp_b*b)+n/(kp_b*b)**2*besselh(n,2,kp_b*b)-n/(kp_b*b)*(besselh(n-1,2,kp_b*b)-n/(kp_b*b)*besselh(n,2,kp_b*b))))-1i*omega*rho_b*vp*kp_b*(besselh(n-1,2,kp_b*b)-n/(kp_b*b)*besselh(n,2,kp_b*b))  
          -2*mu_b*n/b*(kp_b*(besselh(n-1,2,kp_b*b)-n/(kp_b*b)*besselh(n,2,kp_b*b))-1/b*besselh(n,2,kp_b*b))+1i*omega*rho_b*vs*n/b*besselh(n,2,kp_b*b)];

Cn(:,4) = [n/a*besselh(n,2,ks_b*a)
           -ks_b*(besselh(n-1,2,ks_a*b)-n/(ks_b*a)*besselh(n,2,ks_b*a))
           2*mu_b*n/a*(-1/a*(1+n)*besselh(n,2,ks_b*a)+ks_b*besselh(n-1,2,ks_b*a))
           2*mu_b/a*(-1/a*(n**2+n-0.5*(ks_b*a)**2)*besselh(n,2,ks_b*a)+ks_b*besselh(n-1,2,ks_b*a))
           2*mu_b*n/b*(-1/b*besselh(n,2,ks_b*b)+ks_b*(besselh(n-1,2,ks_b*b)-n/(ks_b*b)*besselh(n,2,ks_b*b)))-1i*omega*rho_b*vp*n/b*besselh(n,2,ks_b*b)
           mu_b*(-n**2/b**2*besselh(n,2,ks_b*b)+1/b*ks_b*(besselh(n-1,2,ks_b*b)-n/(ks_b*b)*besselh(n,2,ks_b*b))-ks_b**2*(besselh(n-2,2,ks_b*b)-(n-1)/(ks_b*b)*besselh(n-1,2,ks_b*b)+n/(ks_b*b)**2*besselh(n,2,ks_b*b)-n/(ks_b*b)*(besselh(n-1,2,ks_b*b)-n/(ks_b*b)*besselh(n,2,ks_b*b))))+1i*omega*rho_b*vs*ks_b*(besselh(n-1,2,ks_b*b)-n/(ks_b*b)*besselh(n,2,ks_b*b))];   

Cn(:,5) = [-kp_a*(-besselj(n+1,kp_a*a)+n/(kp_a*a)*besselj(n,kp_a*a))
           n/a*besselj(n,kp_a*a)
           -kp_a**2*((-lambda_a-2*mu_a/(kp_a*a)**2*(n+(kp_a*a)**2-n**2))*besselj(n,kp_a*a)+2*mu_a/(kp_a*a)*besselj(n+1,kp_a*a))
           2*mu_a*n/a*(-kp_a*besselj(n+1,kp_a*a)+(n-1)/a*besselj(n,kp_a*a))
           0
           0];

Cn(:,6) = [-n/a*besselj(n,ks_a*a)
            ks_a*(-besselj(n+1,ks_a*a)+n/(ks_a*a)*besselj(n,ks_a*a))
            -2*mu_a*n/a*(-1/a*(1-n)*besselj(n,ks_a*a)-ks_a*besselj(n+1,ks_a*a))
            -mu_a*2/a*(1/a*(-n**2+n+0.5*(ks_a*a)**2)*besselj(n,ks_a*a)-ks_a*besselj(n+1,ks_a*a))
            0
            0]; 
      
Bn      = [-epsilon(n)*1i**n*kp_b*(-besselj(n+1,kp_b*a)+n/(kp_b*a)*besselj(n,kp_b*a))
           n/a*epsilon(n)*1i**n*besselj(n,kp_b*a)
           -epsilon(n)*1i**n*((-lambda_b*kp_b**2-2*mu_b/a**2*(n-n**2+(kp_b*a)**2))*besselj(n,kp_b*a)+2*mu_b*kp_b/a*besselj(n+1,kp_b*a))
           2*mu_b*epsilon(n)*1i**n*n/a*(-(1-n)/a*besselj(n,kp_b*a)-kp_b*besselj(n+1,kp_b*a))
           0
           0]; 
      
An =  Cn\Bn;
endfunction

function epsilon = epsilon(n)
  if n==0
    epsilon = 1;
  else
    epsilon = 2;
  end
endfunction  