function A=An(dim)
global a b rho w mu vs vp ks kp kp2 ks2 b2


Cn(:,1)=[kp*(-besselh(dim+1,1,kp*a)+(dim/(kp*a))*besselh(dim,1,kp*a))
    -(dim/a)*besselh(dim,1,kp*a)
    (dim*dim-dim-0.5*ks2*b2)*besselh(dim,1,kp*b)+kp*b*besselh(dim+1,1,kp*b)-1i*w*w*rho*(b2/(2*mu))*(-besselh(dim+1,1,kp*b)+(dim/(kp*b))*besselh(dim,1,kp*b))
    -dim*((dim-1)*besselh(dim,1,kp*b)-kp*b*besselh(dim+1,1,kp*b))+1i*w*rho*vs*dim*(b/(2*mu))*besselh(dim,1,kp*b)];

Cn(:,2)=[(dim/a)*besselh(dim,1,ks*a)
    -ks*(-besselh(dim+1,1,ks*a)+(dim/(ks*a))*besselh(dim,1,ks*a))
    dim*((dim-1)*besselh(dim,1,ks*b)-ks*b*besselh(dim+1,1,ks*b))-1i*w*rho*vp*dim*(b/(2*mu))*besselh(dim,1,ks*b)
    -(dim*dim-dim-0.5*ks2*b2)*besselh(dim,1,ks*b)-ks*b*besselh(dim+1,1,ks*b)+1i*w*w*rho*(b2/(2*mu))*(-besselh(dim+1,1,ks*b)+(dim/(ks*b))*besselh(dim,1,ks*b))];

Cn(:,3)=[kp*(besselh(dim-1,2,kp*a)-(dim/(kp*a))*besselh(dim,2,kp*a))
    -(dim/a)*besselh(dim,2,kp*a)
    (-((ks2-2*kp2)*besselh(dim,2,kp*b)-2*kp2*(besselh(dim-2,2,kp*b)-((dim-1)/(kp*b))*besselh(dim-1,2,kp*b)+(dim/(kp2*b2))*besselh(dim,2,kp*b)-(dim/(kp*b))*(besselh(dim-1,2,kp*b)-(dim/(kp*b))*besselh(dim,2,kp*b))))-1i*w*w*rho*(besselh(dim-1,2,kp*b)-(dim/(kp*b))*besselh(dim,2,kp*b)))*b2/2
    -dim*b*(kp*(besselh(dim-1,2,kp*b)-(dim/(kp*b))*besselh(dim,2,kp*b))-(1/b)*besselh(dim,2,kp*b))+1i*w*rho*vs*dim*(b/(2*mu))*besselh(dim,2,kp*b)];

Cn(:,4)=[(dim/a)*besselh(dim,2,ks*a)
    -ks*(besselh(dim-1,2,ks*a)-(dim/(ks*a))*besselh(dim,2,ks*a))
    dim*b*(ks*(besselh(dim-1,2,ks*b)-(dim/(ks*b))*besselh(dim,2,ks*b))-(1/b)*besselh(dim,2,ks*b))-1i*w*rho*vp*dim*(b/(2*mu))*besselh(dim,2,ks*b)
    ((-(dim*dim/b2)*besselh(dim,2,ks*b)+(ks/b)*(besselh(dim-1,2,ks*b)-(dim/(ks*b))*besselh(dim,2,ks*b))-ks2*(besselh(dim-2,2,ks*b)-((dim-1)/(ks*b))*besselh(dim-1,2,ks*b)+(dim/(ks2*b2))*besselh(dim,2,ks*b)-(dim/(ks*b))*(besselh(dim-1,2,ks*b)-(dim/(ks*b))*besselh(dim,2,ks*b))))+1i*w*w*(rho/mu)*(besselh(dim-1,2,ks*b)-(dim/(ks*b))*besselh(dim,2,ks*b)))*b2/2];

Bn=[-epsilon(dim)*(1i^dim)*kp*((dim/(kp*a))*besselj(dim,kp*a)-besselj(dim+1,kp*a))
    (dim/a)*epsilon(dim)*(1i^dim)*besselj(dim,kp*a)
    0
    0];

A=Cn\Bn;
endfunction

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