%
disp(' annular_plate  ver 1.4  November 19, 2009 ');
disp(' by Tom Irvine ');
disp(' ');
disp(' This program calculates the natural frequencies ');
disp(' of the bending modes of an annular plate ');
disp(' Reference: Leissa, Vibration of Plates NASA SP-160 ');
%
disp(' ');
disp(' Enter unit type: 1=English 2=metric');
iu=input(' ');
%
if(iu==1)
    disp(' ');
    out_diam=input(' Enter outside diameter (inch) ');
    disp(' ');
    in_diam=input(' Enter inside diameter (inch) ');
    disp(' '); 
    h=input(' Enter thickness (inch) ');
else
    disp(' ');
    disp(' Select diameter unit:');
    disp(' 1=meters  2=cm  3=mm ');
    ilu=input(' ');  
    if(ilu==1)
        disp(' ');
        out_diam=input(' Enter outside diameter (m) ');  
        disp(' ');
        in_diam=input(' Enter inside diameter (m) ');      
    end
    if(ilu==2)
        disp(' ');
        out_diam=input(' Enter outside diameter (cm) ');
        disp(' ');
        in_diam=input(' Enter inside diameter (cm) ');    
        out_diam=out_diam/100;
        in_diam=in_diam/100;    
    end  
    if(ilu==3)
        disp(' ');
        out_diam=input(' Enter outside diameter (mm) ');
        disp(' ');
        in_diam=input(' Enter inside diameter (mm) ');   
        out_diam=out_diam/1000;
        in_diam=in_diam/1000;       
    end   
    disp(' ');
    h=input(' Enter thickness (mm) ');
    h=h/1000;
end
%
disp(' ');
imat=input(' Enter material: \n 1=aluminum  2=steel  3=G10  4=other ');
%
if(imat==1) % Al
    if(iu==1)
        rho=0.1;
        E=1.0e+07;
    else
        rho=2768;
        E=6.891e+010;
    end
    mu=0.3;
end
if(imat==2) % steel
    if(iu==1)
        rho=0.29;
        E=3.0e+007;
    else
       rho=7887;
       E=2.067e+011;
    end
    mu=0.3;
end
if(imat==3) % G10
    if(iu==1)
        rho=0.065;
        E=2.7e+006;
    else
       rho=1799;
       E=1.861e+010;        
    end
    mu=0.12;
end
if(imat ~=1 & imat ~=2 & imat ~=3)
    if(iu==1)
        disp(' ');   
        E = input(' Enter the elastic modulus (lbf/in^2) ');
        disp(' ');
        rho = input(' Enter the mass density (lbm/in^3) '); 
    else
        disp(' ');   
        E = input(' Enter the elastic modulus (N/m^2) ');
        disp(' ');
        rho = input(' Enter the mass density (kg/m^3) ');        
    end
    disp(' ');
    mu=input(' Enter Poisson ratio ');    
end
%
D=E*h^3/(12*(1.-mu^2));
%
smass=rho*(out_diam^2-in_diam^2)*h*pi/4;
disp(' ');
if(iu==1)
    out1=sprintf(' Structural mass = %8.4g lbm',smass);
    disp(out1);
else
    out1=sprintf(' Structural mass = %8.4g kg',smass);
    disp(out1);
end
disp(' ');
disp(' Add uniformly distributed non-structural mass ?');
disp(' 1=yes  2=no ')
imass=input(' ');
%
NSM=0;
if(imass==1)
  disp(' ');
  if(iu==1) 
    disp(' Enter added mass(lbm) ');
  else
    disp(' Enter added mass(kg) ');
  end
  NSM=input(' ');
end
%
mass=NSM+smass;
M=mass;
%
if(iu==1)
    disp(' ');
    out1=sprintf(' Total mass = %8.4g lbm ',M);
    disp(out1);
    %
    disp(' ');
    out1=sprintf(' Plate Stiffness Factor D = %8.4g lbf in',D);
    disp(out1); 
    %
    mass=mass/386;
    M=M/386;
    rho=rho/386;   % mass per volume
else
    disp(' ');
    out1=sprintf(' Total mass = %8.4g kg ',M);
    disp(out1);
    %
    disp(' ');
    out1=sprintf(' Plate Stiffness Factor D = %8.4g N m',D);
    disp(out1);     
end
%
disp(' ');
disp(' Enter boundary condition:')
disp('    outside  inside ');
disp(' 1=  fixed   fixed ');
disp(' 2=  fixed   ss ');
disp(' 3=  fixed   free ');
disp(' 4=  ss      fixed ');
disp(' 5=  ss      ss ');
disp(' 6=  ss      free ');
disp(' 7=  free    fixed ');
disp(' 8=  free    ss ');
disp(' 9=  free    free ');
bc=input(' ');
iflag=0;
%
DD=sqrt(D/(rho*h));
% Bessel Functions of the First Kind
% JN =besselj(n,x);
% Modified Bessel Functions of the First Kind
% IN =besseli(n,x);
% Bessel Functions of the Second Kind
% YN=bessely(n,x);
% Modified Bessel Functions of the Second Kind
% KN=besselk(n,x);
%
    
b = in_diam/2;
a = out_diam/2;
alpha=b/a;
%
dx=0.005;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear J0A;
clear J1A;
clear I0A;
clear I1A;
clear Y0A;
clear Y1A;
clear K0A;
clear K1A;
%
disp(' ')
disp(' Calculating Bessel functions... ');
%
N=4000;
for(i=1:N)
%
    x=i*dx;
    xa=alpha*x;
%
    J0(i)=besselj(0,x);
    J1(i)=besselj(1,x);
%
    I0(i)=besseli(0,x);
    I1(i)=besseli(1,x);
%
    Y0(i)=bessely(0,x);
    Y1(i)=bessely(1,x);
%
    K0(i)=besselk(0,x);
    K1(i)=besselk(1,x);    
%    
    J0A(i)=besselj(0,xa);
    J1A(i)=besselj(1,xa);
%
    I0A(i)=besseli(0,xa);
    I1A(i)=besseli(1,xa);
%
    Y0A(i)=bessely(0,xa);
    Y1A(i)=bessely(1,xa);
%
    K0A(i)=besselk(0,xa);
    K1A(i)=besselk(1,xa);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
if(bc==1)  % fixed-fixed
    iflag=1;
    clear L2;
    n=0;  %bc=1
    k=0;
    ijk=1;
    clear m;
    clear d;
    clear nn;
    clear kk;
    for(i=1:N)
        x=i*dx;
%
        m(1,1)=J0(i);
        m(1,2)=Y0(i);
        m(1,3)=I0(i);
        m(1,4)=K0(i);
%  
        m(2,1)=J1(i);
        m(2,2)=Y1(i);
        m(2,3)=-I1(i);
        m(2,4)=K1(i);
%
        m(3,1)=J0A(i);
        m(3,2)=Y0A(i);
        m(3,3)=I0A(i);
        m(3,4)=K0A(i);
%
        m(4,1)=J1A(i);
        m(4,2)=Y1A(i);
        m(4,3)=-I1A(i);
        m(4,4)=K1A(i);
%    
        dd=det(m);
        d(i)=dd;
        if(i>1)
            if((dd*dbefore)<=0)
                rt=(x+(x-dx))/2;
                out1=sprintf('lambda^2 = %8.5g n=%d',rt^2,n);
%               disp(out1);
                L2(ijk)=rt^2;
                nn(ijk)=n;
                kk(ijk)=k;
                k=k+1;    
                ijk=ijk+1;
            end
        end
        dbefore=dd;
% 
    end
%
    n=1;  %bc=1
    clear m;
    clear d;
    k=0;
    for(i=1:N)
        x=i*dx;
%
        m(1,1)=J1(i);
        m(1,2)=Y1(i);
        m(1,3)=I1(i);
        m(1,4)=K1(i);
%  
        m(2,1)=J0(i);
        m(2,2)=Y0(i);
        m(2,3)=I0(i);
        m(2,4)=-K0(i);
%
        m(3,1)=J1A(i);
        m(3,2)=Y1A(i);
        m(3,3)=I1A(i);
        m(3,4)=K1A(i);
%
        m(4,1)=J0A(i);
        m(4,2)=Y0A(i);
        m(4,3)=I0A(i);
        m(4,4)=-K0A(i);
%    
        dd=det(m);
        d(i)=dd;
        if(i>1)
            if((dd*dbefore)<=0)
                rt=(x+(x-dx))/2;
                out1=sprintf('lambda^2 = %8.5g n=%d',rt^2,n);
                disp(out1);
                L2(ijk)=rt^2;
                nn(ijk)=n;
                kk(ijk)=k;
                k=k+1;     
                ijk=ijk+1;   
            end
        end
        dbefore=dd;
% 
    end
%
    n=2;  %bc=1
    k=0;
    clear m;
    clear d;
    for(i=1:N)
        x=i*dx;
        xa=alpha*x; 
%
        m(1,1)=J0(i);
        m(1,2)=Y0(i);
        m(1,3)=-I0(i)+(4/x)*I1(i);
        m(1,4)=-K0(i)-(4/x)*K1(i);
%  
        m(2,1)=J1(i);
        m(2,2)=Y1(i);
        m(2,3)=I1(i);
        m(2,4)=-K1(i);
%
        m(3,1)=J0A(i);
        m(3,2)=Y0A(i);
        m(3,3)=-I0A(i)+(4/xa)*I1A(i);
        m(3,4)=-K0A(i)-(4/xa)*K1A(i);
%
        m(4,1)=J1A(i);
        m(4,2)=Y1A(i);
        m(4,3)=I1A(i);
        m(4,4)=-K1A(i);
%    
        dd=det(m);
        d(i)=dd;
        if(i>1)
            if((dd*dbefore)<=0)
                rt=(x+(x-dx))/2;
                out1=sprintf('lambda^2 = %8.5g n=%d',rt^2,n);
%               disp(out1);
                L2(ijk)=rt^2;
                nn(ijk)=n; 
                kk(ijk)=k;
                k=k+1;
                ijk=ijk+1;    
            end
        end
        dbefore=dd;
% 
    end
end % end BC 1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
if(bc==2)  % fixed-ss
    iflag=1;
    clear L2;
    n=0;  %bc=2
    k=0;
    ijk=1;
    clear m;
    clear d;
    clear nn;
    clear kk;
    for(i=1:N)
        x=i*dx;
%
        m(1,1)=J0(i);
        m(1,2)=Y0(i);
        m(1,3)=I0(i);
        m(1,4)=K0(i);
%  
        m(2,1)=J1(i);
        m(2,2)=Y1(i);
        m(2,3)=-I1(i);
        m(2,4)=K1(i);
%
        A=-2*xa/(1-mu);
%
        m(3,1)=J0A(i);
        m(3,2)=Y0A(i);
        m(3,3)=I0A(i);
        m(3,4)=K0A(i);
%  
        m(4,1)=-J1A(i);
        m(4,2)=-Y1A(i);
        m(4,3)= I1A(i)+A*I0A(i);
        m(4,4)=-K1A(i)+A*K0A(i);
%    
        dd=det(m);
        d(i)=dd;
        if(i>1)
            if((dd*dbefore)<=0)
                rt=(x+(x-dx))/2;
                out1=sprintf('lambda^2 = %8.5g n=%d',rt^2,n);
%               disp(out1);
                L2(ijk)=rt^2;
                nn(ijk)=n;
                kk(ijk)=k;
                k=k+1;    
                ijk=ijk+1;
            end
        end
        dbefore=dd;
% 
    end
%
    n=1;  %bc=2
    clear m;
    clear d;
    k=0;
    for(i=1:N)
        x=i*dx;
%
        m(1,1)=J1(i);
        m(1,2)=Y1(i);
        m(1,3)=I1(i);
        m(1,4)=K1(i);
%  
        m(2,1)=J0(i);
        m(2,2)=Y0(i);
        m(2,3)=I0(i);
        m(2,4)=-K0(i);
%
        E=2*xa/(1-mu);   
        m(3,1)=J0A(i);
        m(3,2)=Y0A(i);
        m(3,3)= I0A(i)-E*I1A(i);
        m(3,4)=-K0A(i)-E*K1A(i);
%
        m(4,1)=J1A(i);
        m(4,2)=Y1A(i);
        m(4,3)=I1A(i);
        m(4,4)=K1A(i); 
%    
        dd=det(m);
        d(i)=dd;
        if(i>1)
            if((dd*dbefore)<=0)
                rt=(x+(x-dx))/2;
                out1=sprintf('lambda^2 = %8.5g n=%d',rt^2,n);
                disp(out1);
                L2(ijk)=rt^2;
                nn(ijk)=n;
                kk(ijk)=k;
                k=k+1;     
                ijk=ijk+1;   
            end
        end
        dbefore=dd;
% 
    end
%
    n=2;  %bc=2
    k=0;
    clear m;
    clear d;
    for(i=1:N)
        x=i*dx;
        xa=alpha*x; 
%
        m(1,1)=J0(i);
        m(1,2)=Y0(i);
        m(1,3)=-I0(i)+(4/x)*I1(i);
        m(1,4)=-K0(i)-(4/x)*K1(i);
%  
        m(2,1)=J1(i);
        m(2,2)=Y1(i);
        m(2,3)=I1(i);
        m(2,4)=-K1(i);
% 
        E=(5-mu)/(1-mu);
        F=(2*xa)/(1-mu);
        G=(4/xa)*(3-mu)/(1-mu);
%
        m(3,1)=J1A(i);
        m(3,2)=Y1A(i);
        m(3,3)= E*I1A(i)-F*I0A(i);
        m(3,4)=-E*K1A(i)-F*K0A(i);
%
        m(4,1)=J0A(i);
        m(4,2)=Y0A(i);
        m(4,3)= G*I1A(i)-E*I0A(i);
        m(4,4)=-G*K1A(i)-E*K0A(i);
%    
        dd=det(m);
        d(i)=dd;
        if(i>1)
            if((dd*dbefore)<=0)
                rt=(x+(x-dx))/2;
                out1=sprintf('lambda^2 = %8.5g n=%d',rt^2,n);
%               disp(out1);
                L2(ijk)=rt^2;
                nn(ijk)=n; 
                kk(ijk)=k;
                k=k+1;
                ijk=ijk+1;    
            end
        end
        dbefore=dd;
% 
    end
end % end BC 2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
if(bc==3)  % fixed-free
    iflag=1;
    clear L2;
    n=0; %bc=3
    k=0;
    ijk=1;
    clear m;
    clear d;
    clear nn;
    clear kk;
    for(i=1:N)
        x=i*dx;
%
        m(1,1)=J0(i);
        m(1,2)=Y0(i);
        m(1,3)=I0(i);
        m(1,4)=K0(i);
%  
        m(2,1)=J1(i);
        m(2,2)=Y1(i);
        m(2,3)=-I1(i);
        m(2,4)=K1(i);
%
        B=2*(1-mu)/(x*alpha);

        m(3,1)=J0A(i);
        m(3,2)=Y0A(i);
        m(3,3)=-I0A(i)+B*I1A(i);
        m(3,4)=-K0A(i)-B*K1A(i);
%
        m(4,1)=J1A(i);
        m(4,2)=Y1A(i);
        m(4,3)=I1A(i);
        m(4,4)=-K1A(i);
%    
        dd=det(m);
        d(i)=dd;
        q(i)=x^2;
        if(i>1)
            if((dd*dbefore)<=0)
                rt=(x+(x-dx))/2;
                out1=sprintf('lambda^2 = %8.5g n=%d',rt^2,n);
%               disp(out1);
                L2(ijk)=rt^2;
                nn(ijk)=n;
                kk(ijk)=k;
                k=k+1;    
                ijk=ijk+1;
            end
        end
        dbefore=dd;
% 
    end
%
    n=1; %bc=3
    clear m;
    clear d;
    k=0;
    for(i=1:N)
        x=i*dx;
        xa=x*alpha;
%
        m(1,1)=J1(i);
        m(1,2)=Y1(i);
        m(1,3)=I1(i);
        m(1,4)=K1(i);
%  
        m(2,1)=J0(i);
        m(2,2)=Y0(i);
        m(2,3)=I0(i);
        m(2,4)=-K0(i);
%
        B=8*(1-mu)/(xa^3);
%
        m(3,1)=J0A(i);
        m(3,2)=Y0A(i);
        m(3,3)= (-1+(xa/2)*B)*I0A(i)-B*I1A(i);
        m(3,4)=-(-1+(xa/2)*B)*K0A(i)-B*K1A(i);
%
        m(4,1)=J1A(i);
        m(4,2)=Y1A(i);
        m(4,3)= (xa^2/4)*B*I0A(i)-(1+(xa/2)*B)*I1A(i);
        m(4,4)=-(xa^2/4)*B*K0A(i)-(1+(xa/2)*B)*K1A(i);
%    
        dd=det(m);
        dq(i)=dd;
        qq(i)=x^2;
        if(i>1)
            if((dd*dbefore)<=0)
                rt=(x+(x-dx))/2;
                out1=sprintf('lambda^2 = %8.5g n=%d',rt^2,n);
 %               disp(out1);
                L2(ijk)=rt^2;
                nn(ijk)=n;
                kk(ijk)=k;
                k=k+1;     
                ijk=ijk+1;   
            end
        end
        dbefore=dd;
% 
    end    
%
    n=2;  %bc=3
    k=0;
    clear m;
    clear d;
    for(i=1:N)
        x=i*dx;
        xa=alpha*x; 
%
        m(1,1)=J0(i);
        m(1,2)=Y0(i);
        m(1,3)=-I0(i)+(4/x)*I1(i);
        m(1,4)=-K0(i)-(4/x)*K1(i);
%  
        m(2,1)=J1(i);
        m(2,2)=Y1(i);
        m(2,3)=I1(i);
        m(2,4)=-K1(i);
%
        mu3=3+mu;
        muq=(1-mu)^2;
        WA=(xa/4)-(mu3/(2*xa)); 
%
        num=48*(1-mu)*xa;
        CSDSden=12*muq-xa^4;
        CS=num/CSDSden;
 %      
        num=12*(1-mu)*(7+mu+xa^2)-xa^4;
        DS=num/CSDSden;
%      
        AS=1-WA*CS;
%      
        BS=((xa/4)+(mu3/(2*xa)))-WA*DS;
%
        m(3,1)=J0A(i);
        m(3,2)=Y0A(i);
        m(3,3)=AS*I0A(i)-BS*I1A(i);
        m(3,4)=AS*K0A(i)+BS*K1A(i);
%
        m(4,1)=J1A(i);
        m(4,2)=Y1A(i);
        m(4,3)=CS*I0A(i)-DS*I1A(i);
        m(4,4)=CS*K0A(i)+DS*K1A(i);
%        
        dd=det(m);
        d(i)=dd;
        if(i>1)
            if((dd*dbefore)<=0)
%   
%  slope check
%
                slope=(dd-dbefore)/dx;
                if(abs(slope)<1000)                
                    rt=(x+(x-dx))/2;
                    out1=sprintf('lambda^2 = %8.5g n=%d',rt^2,n);
%                   disp(out1);
                    L2(ijk)=rt^2;
                    nn(ijk)=n; 
                    kk(ijk)=k;
                    k=k+1;
                    ijk=ijk+1;    
                end
            end
        end
        dbefore=dd;
% 
    end  
%
end % end BC 3
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(bc==4)  % simply supported - fixed
    iflag=1;
    clear L2;
    n=0;  %bc=4
    k=0;
    ijk=1;
    clear m;
    clear d;
    clear nn;
    clear kk;
    for(i=1:N)
        x=i*dx;
        A=2*x/(1-mu);
%
        m(1,1)=J0(i);
        m(1,2)=Y0(i);
        m(1,3)=I0(i);
        m(1,4)=K0(i);
%  
        m(2,1)=J1(i);
        m(2,2)=Y1(i);
        m(2,3)= A*I0(i)-I1(i);
        m(2,4)=-A*K0(i)-K1(i);
%
        m(3,1)=J0A(i);
        m(3,2)=Y0A(i);
        m(3,3)=I0A(i);
        m(3,4)=K0A(i);
%
        m(4,1)=J1A(i);
        m(4,2)=Y1A(i);
        m(4,3)=-I1A(i);
        m(4,4)=K1A(i);
%    
        dd=det(m);
        d(i)=dd;
        if(i>1)
            if((dd*dbefore)<=0)
                rt=(x+(x-dx))/2;
                out1=sprintf('lambda^2 = %8.5g n=%d',rt^2,n);
%               disp(out1);
                L2(ijk)=rt^2;
                nn(ijk)=n;
                kk(ijk)=k;
                k=k+1;    
                ijk=ijk+1;
            end
        end
        dbefore=dd;
% 
    end
%
    n=1;  %bc=4
    clear m;
    clear d;
    k=0;
    for(i=1:N)
        x=i*dx;
        A=2*x/(1-mu);        
%
        m(1,1)=J1(i);
        m(1,2)=Y1(i);
        m(1,3)=I1(i);
        m(1,4)=K1(i);
%  
        m(2,1)=J0(i);
        m(2,2)=Y0(i);
        m(2,3)= I0(i)-A*I1(i);
        m(2,4)=-K0(i)-A*K1(i);
%
        m(3,1)=J1A(i);
        m(3,2)=Y1A(i);
        m(3,3)=I1A(i);
        m(3,4)=K1A(i);
%
        m(4,1)=J0A(i);
        m(4,2)=Y0A(i);
        m(4,3)=I0A(i);
        m(4,4)=-K0A(i);
%    
        dd=det(m);
        d(i)=dd;
        if(i>1)
            if((dd*dbefore)<=0)
                rt=(x+(x-dx))/2;
                out1=sprintf('lambda^2 = %8.5g n=%d',rt^2,n);
                disp(out1);
                L2(ijk)=rt^2;
                nn(ijk)=n;
                kk(ijk)=k;
                k=k+1;     
                ijk=ijk+1;   
            end
        end
        dbefore=dd;
% 
    end
%
    n=2;  %bc=4
    k=0;
    clear m;
    clear d;
    for(i=1:N)
        x=i*dx;
        xa=alpha*x; 
        A=(5-mu)/(1-mu);
        B=(3-mu)/(1-mu);
        C=2*x/(1-mu);       
%
        m(1,1)=J1(i);
        m(1,2)=Y1(i);
        m(1,3)= A*I1(i)-C*I0(i);
        m(1,4)=-A*K1(i)-C*K0(i);
%
        m(2,1)=J0(i);
        m(2,2)=Y0(i);
        m(2,3)= (4/x)*B*I1(i)-A*I0(i);
        m(2,4)=-(4/x)*B*K1(i)-A*K0(i);
%
        m(3,1)=J0A(i);
        m(3,2)=Y0A(i);
        m(3,3)=-I0A(i)+(4/xa)*I1A(i);
        m(3,4)=-K0A(i)-(4/xa)*K1A(i);
%
        m(4,1)=J1A(i);
        m(4,2)=Y1A(i);
        m(4,3)=I1A(i);
        m(4,4)=-K1A(i);
%    
        dd=det(m);
        d(i)=dd;
        if(i>1)
            if((dd*dbefore)<=0)
                rt=(x+(x-dx))/2;
                out1=sprintf('lambda^2 = %8.5g n=%d',rt^2,n);
%               disp(out1);
                L2(ijk)=rt^2;
                nn(ijk)=n; 
                kk(ijk)=k;
                k=k+1;
                ijk=ijk+1;    
            end
        end
        dbefore=dd;
% 
    end
end % end BC 4
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(bc==5) % simply-suppored both edges
    iflag=1;
    clear L2;
    n=0;  %bc=5
    k=0;
    ijk=1;
    clear m;
    clear d;
    clear nn;
    clear kk;
    for(i=1:N)
        x=i*dx;
%
        A=-2*x/(1-mu);
%
        m(1,1)=J0(i);
        m(1,2)=Y0(i);
        m(1,3)=I0(i);
        m(1,4)=K0(i);
%  
        m(2,1)=-J1(i);
        m(2,2)=-Y1(i);
        m(2,3)= I1(i)+A*I0(i);
        m(2,4)=-K1(i)+A*K0(i);
%
        A=-2*xa/(1-mu);
%
        m(3,1)=J0A(i);
        m(3,2)=Y0A(i);
        m(3,3)=I0A(i);
        m(3,4)=K0A(i);
%  
        m(4,1)=-J1A(i);
        m(4,2)=-Y1A(i);
        m(4,3)= I1A(i)+A*I0A(i);
        m(4,4)=-K1A(i)+A*K0A(i);
%    
        dd=det(m);
        d(i)=dd;
        if(i>1)
            if((dd*dbefore)<=0)
                rt=(x+(x-dx))/2;
                out1=sprintf('lambda^2 = %8.5g n=%d',rt^2,n);
%               disp(out1);
                L2(ijk)=rt^2;
                nn(ijk)=n;
                kk(ijk)=k;
                k=k+1;    
                ijk=ijk+1;
            end
        end
        dbefore=dd;
% 
    end
%
    n=1;  %bc=5
    clear m;
    clear d;
    k=0;
    for(i=1:N)
%%
        x=i*dx;
        xa=x*alpha;
        E=2*x/(1-mu);        
%  
        m(1,1)=J0(i);
        m(1,2)=Y0(i);
        m(1,3)= I0(i)-E*I1(i);
        m(1,4)=-K0(i)-E*K1(i);
%
        m(2,1)=J1(i);
        m(2,2)=Y1(i);
        m(2,3)=I1(i);
        m(2,4)=K1(i);        
%
        E=2*xa/(1-mu);   
        m(3,1)=J0A(i);
        m(3,2)=Y0A(i);
        m(3,3)= I0A(i)-E*I1A(i);
        m(3,4)=-K0A(i)-E*K1A(i);
%
        m(4,1)=J1A(i);
        m(4,2)=Y1A(i);
        m(4,3)=I1A(i);
        m(4,4)=K1A(i); 
%    
%%
        dd=det(m);
        d(i)=dd;
        if(i>1)
            if((dd*dbefore)<=0)
                rt=(x+(x-dx))/2;
                out1=sprintf('lambda^2 = %8.5g n=%d',rt^2,n);
%                disp(out1);
                L2(ijk)=rt^2;
                nn(ijk)=n;
                kk(ijk)=k;
                k=k+1;     
                ijk=ijk+1;  
                if(k==5)
                    break;
                end
            end
        end
        dbefore=dd;
% 
    end
%
    n=2;  %bc=5
    k=0;
    clear m;
    clear d;
    for(i=1:N)
        x=i*dx;
        xa=alpha*x;  
%
        A=-((3+mu)/(2*xa))+(xa/4);
        B= ((3+mu)/(2*xa))+(xa/4);
%       
        num=48*(1-mu)*xa;
        den=12*muq-xa^4;
        C=num/den;
        num=12*(1-mu)*(7+mu+xa^2)-xa^4;
        D=num/den;
%
        E=(5-mu)/(1-mu);
        F=(2*x)/(1-mu);
        G=(4/x)*(3-mu)/(1-mu);
        AS=1-A*C;
        BS=B-A*D;
%
        m(1,1)=J1(i);
        m(1,2)=Y1(i);
        m(1,3)= E*I1(i)-F*I0(i);
        m(1,4)=-E*K1(i)-F*K0(i);
%
        m(2,1)=J0(i);
        m(2,2)=Y0(i);
        m(2,3)= G*I1(i)-E*I0(i);
        m(2,4)=-G*K1(i)-E*K0(i);
%
        F=(2*xa)/(1-mu);
        G=(4/xa)*(3-mu)/(1-mu);
%
        m(3,1)=J1A(i);
        m(3,2)=Y1A(i);
        m(3,3)= E*I1A(i)-F*I0A(i);
        m(3,4)=-E*K1A(i)-F*K0A(i);
%
        m(4,1)=J0A(i);
        m(4,2)=Y0A(i);
        m(4,3)= G*I1A(i)-E*I0A(i);
        m(4,4)=-G*K1A(i)-E*K0A(i);
%    
        dd=det(m);
        d(i)=dd;
        d72(i)=dd;
        q72(i)=x^2;      
        if(i>1)
            if((dd*dbefore)<=0)
                slope=(dd-dbefore)/dx;
                if(abs(slope)<1000 && x^2 > 20.)  
                    rt=(x+(x-dx))/2;
                    out1=sprintf('lambda^2 = %8.5g n=%d',rt^2,n);
%                   disp(out1);
                    L2(ijk)=rt^2;
                    nn(ijk)=n; 
                    kk(ijk)=k;
                    k=k+1;
                    ijk=ijk+1;
                end
            end
        end
        dbefore=dd;
% 
    end
end % end BC 5
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(bc==6) % simply-suppored free 
    iflag=1;
    clear L2;
    n=0;  %bc=6
    k=0;
    ijk=1;
    clear m;
    clear d;
    clear nn;
    clear kk;
    for(i=1:N)
        x=i*dx;
        xa=x*alpha;
        A=-2*x/(1-mu);
%
        m(1,1)=J0(i);
        m(1,2)=Y0(i);
        m(1,3)=I0(i);
        m(1,4)=K0(i);
%  
        m(2,1)=-J1(i);
        m(2,2)=-Y1(i);
        m(2,3)= I1(i)+A*I0(i);
        m(2,4)=-K1(i)+A*K0(i);
%
        B=2*(1-mu)/(x*alpha);
%
        m(3,1)=J0A(i);
        m(3,2)=Y0A(i);
        m(3,3)=-I0A(i)+B*I1A(i);
        m(3,4)=-K0A(i)-B*K1A(i);
%
        m(4,1)=J1A(i);
        m(4,2)=Y1A(i);
        m(4,3)=I1A(i);
        m(4,4)=-K1A(i);
%    
        dd=det(m);
        d(i)=dd;
        if(i>1)
            if((dd*dbefore)<=0)
                rt=(x+(x-dx))/2;
                out1=sprintf('lambda^2 = %8.5g n=%d',rt^2,n);
%               disp(out1);
                L2(ijk)=rt^2;
                nn(ijk)=n;
                kk(ijk)=k;
                k=k+1;    
                ijk=ijk+1;
            end
        end
        dbefore=dd;
% 
    end
%
    n=1;  %bc=6
    clear m;
    clear d;
    k=0;
    for(i=1:N)
%% 
        x=i*dx;
        xa=x*alpha;
        E=2*x/(1-mu);        
%  
        m(1,1)=J0(i);
        m(1,2)=Y0(i);
        m(1,3)= I0(i)-E*I1(i);
        m(1,4)=-K0(i)-E*K1(i);
%
        m(2,1)=J1(i);
        m(2,2)=Y1(i);
        m(2,3)=I1(i);
        m(2,4)=K1(i);  
%    
        B=8*(1-mu)/(xa^3);
%
        m(3,1)=J0A(i);
        m(3,2)=Y0A(i);
        m(3,3)= (-1+(xa/2)*B)*I0A(i)-B*I1A(i);
        m(3,4)=-(-1+(xa/2)*B)*K0A(i)-B*K1A(i);
%
        m(4,1)=J1A(i);
        m(4,2)=Y1A(i);
        m(4,3)= (xa^2/4)*B*I0A(i)-(1+(xa/2)*B)*I1A(i);
        m(4,4)=-(xa^2/4)*B*K0A(i)-(1+(xa/2)*B)*K1A(i);
%%
        dd=det(m);
        d(i)=dd;
        if(i>1)
            if((dd*dbefore)<=0)
                rt=(x+(x-dx))/2;
                out1=sprintf('lambda^2 = %8.5g n=%d',rt^2,n);
%                disp(out1);
                L2(ijk)=rt^2;
                nn(ijk)=n;
                kk(ijk)=k;
                k=k+1;     
                ijk=ijk+1;  
                if(k==5)
                    break;
                end
            end
        end
        dbefore=dd;
% 
    end
%
    n=2;  %bc=6 
    k=0;
    clear m;
    clear d;
    for(i=1:N)
        x=i*dx;
        xa=alpha*x;
%
        A=-((3+mu)/(2*xa))+(xa/4);
        B= ((3+mu)/(2*xa))+(xa/4);
%       
        num=48*(1-mu)*xa;
        den=12*muq-xa^4;
        C=num/den;
        num=12*(1-mu)*(7+mu+xa^2)-xa^4;
        D=num/den;
%
        E=(5-mu)/(1-mu);
        F=(2*x)/(1-mu);
        G=(4/x)*(3-mu)/(1-mu);
        AS=1-A*C;
        BS=B-A*D;
%
        m(1,1)=J1(i);
        m(1,2)=Y1(i);
        m(1,3)= E*I1(i)-F*I0(i);
        m(1,4)=-E*K1(i)-F*K0(i);
%
        m(2,1)=J0(i);
        m(2,2)=Y0(i);
        m(2,3)= G*I1(i)-E*I0(i);
        m(2,4)=-G*K1(i)-E*K0(i);
% 
        mu3=3+mu;
        muq=(1-mu)^2;
        W=(x/4)-(mu3/(2*x));
        WA=(xa/4)-(mu3/(2*xa)); 
%        
        num=48*(1-mu)*x;
        CDden=12*muq-x^4;
        C=num/CDden;
%       
        num=12*(1-mu)*(7+mu+x^2)-x^4;
        D=num/CDden;
 %       
        A=1-W*C;
 %      
        B=((x/4)+(mu3/(2*x)))-W*D;
 %      
        num=48*(1-mu)*xa;
        CSDSden=12*muq-xa^4;
        CS=num/CSDSden;
 %      
        num=12*(1-mu)*(7+mu+xa^2)-xa^4;
        DS=num/CSDSden;
 %      
        AS=1-WA*CS;
%      
        BS=((xa/4)+(mu3/(2*xa)))-WA*DS;
%
        m(3,1)=J0A(i);
        m(3,2)=Y0A(i);
        m(3,3)=AS*I0A(i)-BS*I1A(i);
        m(3,4)=AS*K0A(i)+BS*K1A(i);
%
        m(4,1)=J1A(i);
        m(4,2)=Y1A(i);
        m(4,3)=CS*I0A(i)-DS*I1A(i);
        m(4,4)=CS*K0A(i)+DS*K1A(i);
%      
        dd=det(m);
        d(i)=dd;
        d72(i)=dd;
        q72(i)=x^2;      
        if(i>1)
            if((dd*dbefore)<=0)
                slope=(dd-dbefore)/dx;
                if(abs(slope)<1000)  
                    rt=(x+(x-dx))/2;
                    out1=sprintf('lambda^2 = %8.5g n=%d',rt^2,n);
%                   disp(out1);
                    L2(ijk)=rt^2;
                    nn(ijk)=n; 
                    kk(ijk)=k;
                    k=k+1;
                    ijk=ijk+1;
                end
            end
        end
        dbefore=dd;
% 
    end
end % end BC 6
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(bc==7)  % free-fixed
    iflag=1;
    clear L2;
    n=0;  %bc=7
    k=0;
    ijk=1;
    clear m;
    clear d;
    clear nn;
    clear kk;
    for(i=1:N)
        x=i*dx;
% 
        A=2*(1-mu)/x;
%
        m(1,1)=J0(i);
        m(1,2)=Y0(i);
        m(1,3)=-I0(i)+A*I1(i);
        m(1,4)=-K0(i)-A*K1(i);
%  
        m(2,1)=J1(i);
        m(2,2)=Y1(i);
        m(2,3)=I1(i);
        m(2,4)=-K1(i);
%
        m(3,1)=J0A(i);
        m(3,2)=Y0A(i);
        m(3,3)=I0A(i);
        m(3,4)=K0A(i);
%
        m(4,1)=J1A(i);
        m(4,2)=Y1A(i);
        m(4,3)=-I1A(i);
        m(4,4)=K1A(i);
%    
        dd=det(m);
        d(i)=dd;
        if(i>1)
            if((dd*dbefore)<=0)
                rt=(x+(x-dx))/2;
                out1=sprintf('lambda^2 = %8.5g n=%d',rt^2,n);
%               disp(out1);
                L2(ijk)=rt^2;
                nn(ijk)=n;
                kk(ijk)=k;
                k=k+1;    
                ijk=ijk+1;
            end
        end
        dbefore=dd;
% 
    end
%
    n=1;  %bc=7
    clear m;
    clear d;
    k=0;
    for(i=1:N)
        x=i*dx;
%
        A=8*(1-mu)/x^3;
%
        m(1,1)=J0(i);
        m(1,2)=Y0(i);
        m(1,3)= (-1+(x/2)*A)*I0(i)-A*I1(i);
        m(1,4)=-(-1+(x/2)*A)*K0(i)-A*K1(i);
%  
        m(2,1)=J1(i);
        m(2,2)=Y1(i);
        m(2,3)= (x^2/4)*A*I0(i)-(1+(x/2)*A)*I1(i);
        m(2,4)=-(x^2/4)*A*K0(i)-(1+(x/2)*A)*K1(i);
%
        m(3,1)=J1A(i);
        m(3,2)=Y1A(i);
        m(3,3)=I1A(i);
        m(3,4)=K1A(i);
%
        m(4,1)=J0A(i);
        m(4,2)=Y0A(i);
        m(4,3)=I0A(i);
        m(4,4)=-K0A(i);
%    
        dd=det(m);
        d(i)=dd;
        if(i>1)
            if((dd*dbefore)<=0)
                rt=(x+(x-dx))/2;
                out1=sprintf('lambda^2 = %8.5g n=%d',rt^2,n);
                disp(out1);
                L2(ijk)=rt^2;
                nn(ijk)=n;
                kk(ijk)=k;
                k=k+1;     
                ijk=ijk+1;   
            end
        end
        dbefore=dd;
% 
    end
%
    n=2;  %bc=7
    k=0;
    clear m;
    clear d;
    for(i=1:N)
        x=i*dx;
        xa=x*alpha;
        mu3=3+mu;
        muq=(1-mu)^2;
        W=(x/4)-(mu3/(2*x));
        WA=(xa/4)-(mu3/(2*xa)); 
%        
        num=48*(1-mu)*x;
        CDden=12*muq-x^4;
        C=num/CDden;
%       
        num=12*(1-mu)*(7+mu+x^2)-x^4;
        D=num/CDden;
 %       
        A=1-W*C;
 %      
        B=((x/4)+(mu3/(2*x)))-W*D;
 %      
        num=48*(1-mu)*xa;
        CSDSden=12*muq-xa^4;
        CS=num/CSDSden;
 %      
        num=12*(1-mu)*(7+mu+xa^2)-xa^4;
        DS=num/CSDSden;
 %      
        AS=1-WA*CS;
%      
        BS=((xa/4)+(mu3/(2*xa)))-WA*DS;
%
        m(1,1)=J0(i);
        m(1,2)=Y0(i);
        m(1,3)=A*I0(i)-B*I1(i);
        m(1,4)=A*K0(i)+B*K1(i);
%  
        m(2,1)=J1(i);
        m(2,2)=Y1(i);
        m(2,3)=C*I0(i)-D*I1(i);
        m(2,4)=C*K0(i)+D*K1(i);
%
        m(3,1)=J0A(i);
        m(3,2)=Y0A(i);
        m(3,3)=-I0A(i)+(4/xa)*I1A(i);
        m(3,4)=-K0A(i)-(4/xa)*K1A(i);
%
        m(4,1)=J1A(i);
        m(4,2)=Y1A(i);
        m(4,3)=I1A(i);
        m(4,4)=-K1A(i);
%    
        dd=det(m);
        d(i)=dd;
        if(i>1)
            if((dd*dbefore)<=0)
                slope=(dd-dbefore)/dx;
                if(abs(slope)<1000)        
                    rt=(x+(x-dx))/2;
                    out1=sprintf('lambda^2 = %8.5g n=%d',rt^2,n);
%                   disp(out1);
                    L2(ijk)=rt^2;
                    nn(ijk)=n; 
                    kk(ijk)=k;
                    k=k+1;
                    ijk=ijk+1;    
                end
            end
        end
        dbefore=dd;
% 
    end
end % end BC 7
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(bc==8) % free simply-suppored 
    iflag=1;
    clear L2;
    n=0;  %bc=8
    k=0;
    ijk=1;
    clear m;
    clear d;
    clear nn;
    clear kk;
    for(i=1:N)
        x=i*dx;
        xa=x*alpha;
%
        A=2*(1-mu)/x;
%
        m(1,1)=J0(i);
        m(1,2)=Y0(i);
        m(1,3)=-I0(i)+A*I1(i);
        m(1,4)=-K0(i)-A*K1(i);
%  
        m(2,1)=J1(i);
        m(2,2)=Y1(i);
        m(2,3)=I1(i);
        m(2,4)=-K1(i);
%
        A=-2*xa/(1-mu);
%
        m(3,1)=J0A(i);
        m(3,2)=Y0A(i);
        m(3,3)=I0A(i);
        m(3,4)=K0A(i);
%  
        m(4,1)=-J1A(i);
        m(4,2)=-Y1A(i);
        m(4,3)= I1A(i)+A*I0A(i);
        m(4,4)=-K1A(i)+A*K0A(i);
%    
        dd=det(m);
        d(i)=dd;
        if(i>1)
            if((dd*dbefore)<=0)
                rt=(x+(x-dx))/2;
                out1=sprintf('lambda^2 = %8.5g n=%d',rt^2,n);
%               disp(out1);
                L2(ijk)=rt^2;
                nn(ijk)=n;
                kk(ijk)=k;
                k=k+1;    
                ijk=ijk+1;
            end
        end
        dbefore=dd;
% 
    end
%
    n=1;  %bc=8
    clear m;
    clear d;
    k=0;
    for(i=1:N)
%%
        x=i*dx;
        xa=x*alpha;
% 
        A=8*(1-mu)/x^3;
%
        m(1,1)=J0(i);
        m(1,2)=Y0(i);
        m(1,3)= (-1+(x/2)*A)*I0(i)-A*I1(i);
        m(1,4)=-(-1+(x/2)*A)*K0(i)-A*K1(i);
%  
        m(2,1)=J1(i);
        m(2,2)=Y1(i);
        m(2,3)= (x^2/4)*A*I0(i)-(1+(x/2)*A)*I1(i);
        m(2,4)=-(x^2/4)*A*K0(i)-(1+(x/2)*A)*K1(i);
%
        E=2*xa/(1-mu);   
        m(3,1)=J0A(i);
        m(3,2)=Y0A(i);
        m(3,3)= I0A(i)-E*I1A(i);
        m(3,4)=-K0A(i)-E*K1A(i);
%
        m(4,1)=J1A(i);
        m(4,2)=Y1A(i);
        m(4,3)=I1A(i);
        m(4,4)=K1A(i); 
%    
%%
        dd=det(m);
        d(i)=dd;
        if(i>1)
            if((dd*dbefore)<=0)
                rt=(x+(x-dx))/2;
                out1=sprintf('lambda^2 = %8.5g n=%d',rt^2,n);
%                disp(out1);
                L2(ijk)=rt^2;
                nn(ijk)=n;
                kk(ijk)=k;
                k=k+1;     
                ijk=ijk+1;  
                if(k==5)
                    break;
                end
            end
        end
        dbefore=dd;
% 
    end
%
    n=2;  %bc=8
    k=0;
    clear m;
    clear d;
    for(i=1:N)
        x=i*dx;
        xa=alpha*x;  
% 
        mu3=3+mu;
        muq=(1-mu)^2;
        W=(x/4)-(mu3/(2*x));
%        
        num=48*(1-mu)*x;
        CDden=12*muq-x^4;
        C=num/CDden;
%       
        num=12*(1-mu)*(7+mu+x^2)-x^4;
        D=num/CDden;
%       
        A=1-W*C;
%      
        B=((x/4)+(mu3/(2*x)))-W*D;
%
        m(1,1)=J0(i);
        m(1,2)=Y0(i);
        m(1,3)=A*I0(i)-B*I1(i);
        m(1,4)=A*K0(i)+B*K1(i);
%  
        m(2,1)=J1(i);
        m(2,2)=Y1(i);
        m(2,3)=C*I0(i)-D*I1(i);
        m(2,4)=C*K0(i)+D*K1(i);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%
        E=(5-mu)/(1-mu);
        F=(2*xa)/(1-mu);
        G=(4/xa)*(3-mu)/(1-mu);
%
        m(3,1)=J1A(i);
        m(3,2)=Y1A(i);
        m(3,3)= E*I1A(i)-F*I0A(i);
        m(3,4)=-E*K1A(i)-F*K0A(i);
%
        m(4,1)=J0A(i);
        m(4,2)=Y0A(i);
        m(4,3)= G*I1A(i)-E*I0A(i);
        m(4,4)=-G*K1A(i)-E*K0A(i);
%    
        dd=det(m);
        d(i)=dd;
        d72(i)=dd;
        q72(i)=x^2;      
        if(i>1)
            if((dd*dbefore)<=0)
                slope=(dd-dbefore)/dx;
                if(abs(slope)<1000)  
                    rt=(x+(x-dx))/2;
                    out1=sprintf('lambda^2 = %8.5g n=%d',rt^2,n);
%                   disp(out1);
                    L2(ijk)=rt^2;
                    nn(ijk)=n; 
                    kk(ijk)=k;
                    k=k+1;
                    ijk=ijk+1;
                end
            end
        end
        dbefore=dd;
% 
    end
end % end BC 8
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
if(bc==9)  % free-free
    iflag=1;
    clear L2;
    n=0;  %bc=9
    k=1;
    ijk=1;
    clear m;
    clear d;
    clear nn;
    clear kk;
    for(i=1:N)
        x=i*dx;
        A=2*(1-mu)/x;
        B=2*(1-mu)/(x*alpha);
%
        m(1,1)=J0(i);
        m(1,2)=Y0(i);
        m(1,3)=-I0(i)+A*I1(i);
        m(1,4)=-K0(i)-A*K1(i);
%  
        m(2,1)=J1(i);
        m(2,2)=Y1(i);
        m(2,3)=I1(i);
        m(2,4)=-K1(i);
%
        m(3,1)=J0A(i);
        m(3,2)=Y0A(i);
        m(3,3)=-I0A(i)+B*I1A(i);
        m(3,4)=-K0A(i)-B*K1A(i);
%
        m(4,1)=J1A(i);
        m(4,2)=Y1A(i);
        m(4,3)=I1A(i);
        m(4,4)=-K1A(i);
%    
        dd=det(m);
        if(i>1)
            if((dd*dbefore)<=0)
                slope=(dd-dbefore)/dx;
                if(abs(slope)<1000)            
                    rt=(x+(x-dx))/2;
                    out1=sprintf('lambda^2 = %8.5g n=%d',rt^2,n);
 %                  disp(out1);
                    L2(ijk)=rt^2;
                    nn(ijk)=n;
                    kk(ijk)=k;
                    k=k+1;    
                    ijk=ijk+1;
                end
            end
        end
        dbefore=dd;
% 
    end
%
    n=1; %bc=9
    clear m;
    clear d;
    k=1;
    for(i=1:N)
        x=i*dx;
        xa=x*alpha;
        A=8*(1-mu)/x^3;
        B=8*(1-mu)/(xa^3);
%
        m(1,1)=J0(i);
        m(1,2)=Y0(i);
        m(1,3)= (-1+(x/2)*A)*I0(i)-A*I1(i);
        m(1,4)=-(-1+(x/2)*A)*K0(i)-A*K1(i);
%  
        m(2,1)=J1(i);
        m(2,2)=Y1(i);
        m(2,3)= (x^2/4)*A*I0(i)-(1+(x/2)*A)*I1(i);
        m(2,4)=-(x^2/4)*A*K0(i)-(1+(x/2)*A)*K1(i);
%
        m(3,1)=J0A(i);
        m(3,2)=Y0A(i);
        m(3,3)= (-1+(xa/2)*B)*I0A(i)-B*I1A(i);
        m(3,4)=-(-1+(xa/2)*B)*K0A(i)-B*K1A(i);
%
        m(4,1)=J1A(i);
        m(4,2)=Y1A(i);
        m(4,3)= (xa^2/4)*B*I0A(i)-(1+(xa/2)*B)*I1A(i);
        m(4,4)=-(xa^2/4)*B*K0A(i)-(1+(xa/2)*B)*K1A(i);
%    
        dd=det(m);
        d(i)=dd;
        if(i>1)
            if((dd*dbefore)<=0)
                rt=(x+(x-dx))/2;
                out1=sprintf('lambda^2 = %8.5g n=%d',rt^2,n);
 %               disp(out1);
                L2(ijk)=rt^2;
                nn(ijk)=n;
                kk(ijk)=k;
                k=k+1;     
                ijk=ijk+1;   
            end
        end
        dbefore=dd;
% 
    end
%
    n=2; %bc=9
    clear m;
    clear d;
    clear q;
    k=0;
    for(i=1:N)
        x=i*dx;
        xa=x*alpha;
        mu3=3+mu;
        muq=(1-mu)^2;
        W=(x/4)-(mu3/(2*x));
        WA=(xa/4)-(mu3/(2*xa)); 
%        
        num=48*(1-mu)*x;
        CDden=12*muq-x^4;
        C=num/CDden;
%       
        num=12*(1-mu)*(7+mu+x^2)-x^4;
        D=num/CDden;
 %       
        A=1-W*C;
 %      
        B=((x/4)+(mu3/(2*x)))-W*D;
 %      
        num=48*(1-mu)*xa;
        CSDSden=12*muq-xa^4;
        CS=num/CSDSden;
 %      
        num=12*(1-mu)*(7+mu+xa^2)-xa^4;
        DS=num/CSDSden;
 %      
        AS=1-WA*CS;
%      
        BS=((xa/4)+(mu3/(2*xa)))-WA*DS;
%
        m(1,1)=J0(i);
        m(1,2)=Y0(i);
        m(1,3)=A*I0(i)-B*I1(i);
        m(1,4)=A*K0(i)+B*K1(i);
%  
        m(2,1)=J1(i);
        m(2,2)=Y1(i);
        m(2,3)=C*I0(i)-D*I1(i);
        m(2,4)=C*K0(i)+D*K1(i);
%
        m(3,1)=J0A(i);
        m(3,2)=Y0A(i);
        m(3,3)=AS*I0A(i)-BS*I1A(i);
        m(3,4)=AS*K0A(i)+BS*K1A(i);
%
        m(4,1)=J1A(i);
        m(4,2)=Y1A(i);
        m(4,3)=CS*I0A(i)-DS*I1A(i);
        m(4,4)=CS*K0A(i)+DS*K1A(i);
%    
        dd=det(m);
        d(i)=dd;
        q(i)=x^2;
        if(i>1)
            if((dd*dbefore)<=0)
%   
%  slope check
%
                slope=(dd-dbefore)/dx;
                if(abs(slope)<1000)
                    rt=(x+(x-dx))/2;
                    out1=sprintf('lambda^2 = %8.5g n=%d',rt^2,n);
%                    disp(out1);
                    L2(ijk)=rt^2;
                    nn(ijk)=n;
                    kk(ijk)=k;
                    k=k+1;     
                    ijk=ijk+1;   
               end
            end
        end
        dbefore=dd;
% 
    end
%
end % end BC 9
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  output
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
if(iflag==1)
    clear omegan;
    clear fn;
    clear big;
%
    s=size(L2);
    if(s(2)>s(1))
        L2=L2';
    end
%
    s=size(nn);
    if(s(2)>s(1))
        nn=nn';
    end
%
    s=size(kk);
    if(s(2)>s(1))
        kk=kk';
    end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
clear Xout;
clear Yout;
clear Xin;
clear Xout;
clear X1;
clear Y1;
clear X2;
clear Y2;
clear X3;
clear Y3;
clear X4;
clear Y4;
clear X5;
clear Y5;
%
NOP=360;
rc=1;
delta = rc*(1-alpha);
THETA=linspace(0,2*pi,NOP);
RHO=ones(1,NOP)*rc;
[Xout,Yout] = pol2cart(THETA,RHO);
RHO=ones(1,NOP)*rc*alpha;
[Xin,Yin] = pol2cart(THETA,RHO);
RHO=ones(1,NOP)*(rc+rc*alpha)*0.5;
[X1,Y1] = pol2cart(THETA,RHO);
RHO=ones(1,NOP)*(rc*alpha+delta*0.333);
[X2,Y2] = pol2cart(THETA,RHO);
RHO=ones(1,NOP)*(rc*alpha+delta*0.667);
[X3,Y3] = pol2cart(THETA,RHO);
RHO=ones(1,NOP)*(rc*alpha+delta*0.25);
[X4,Y4] = pol2cart(THETA,RHO);
RHO=ones(1,NOP)*(rc*alpha+delta*0.75);
[X5,Y5] = pol2cart(THETA,RHO);
%
D1X(1)=0;
D1X(2)=0;
D1pY(1)=rc;
D1pY(2)=rc*alpha;
D1nY(1)=-rc;
D1nY(2)=-rc*alpha;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    disp(' ');
    disp(' Eigenvalues   n    k');
    for(i=1:max(size(L2)))
        out1=sprintf('%9.4g \t  %d \t%d',L2(i),nn(i),kk(i));
        disp(out1);
    end
%
    big=[L2,nn,kk];
    big=sortrows(big,1);
%
    L2=big(:,1);
    nn=big(:,2);
    kk=big(:,3);
%
    omegan=L2*DD/(a^2);
    fn=omegan/(2*pi);
    disp(' ');
    disp('  Natural ');
    disp(' Frequency ');
    disp('   (Hz)       n    k');
    for(i=1:max(size(L2)))
        out1=sprintf('%9.4g \t  %d \t%d',fn(i),nn(i),kk(i));
        disp(out1);
        if(nn(i)<=2 && kk(i)<=3)
            if(i==1)
                figure(1);
            end
            if(i==2)
                figure(2);
            end  
            if(i==3)
                figure(3);
            end
            if(i==4)
                figure(4);
            end 
            if(i==5)
                figure(5);
            end       
            if(i==6)
                figure(6);
            end  
            if(i==7)
                figure(7);
            end     
            if(i==8)
                figure(8);
            end     
            if(i==9)
                figure(9);
            end  
            if(i==10)
                figure(10);
            end     
            hold on;
            plot(Xout,Yout);
            plot(Xin,Yin);  
            out1=sprintf(' Mode %d  fn=%8.4g Hz  n=%d k=%d',i,fn(i),nn(i),kk(i));
            title(out1);
            if(nn(i)==1)
                plot(D1X,D1pY);
                plot(D1X,D1nY);
            end
            if(nn(i)==2)
                plot(D1X,D1pY);
                plot(D1X,D1nY); 
                plot(D1pY,D1X);
                plot(D1nY,D1X);          
            end    
            if(kk(i)==1)
                plot(X1,Y1);               
            end  
            if(kk(i)==2)
                plot(X2,Y2);
                plot(X3,Y3);    
            end  
            if(kk(i)==3)
                plot(X1,Y1);
                plot(X4,Y4);  
                plot(X5,Y5);        
            end  
            axis square;
            hold off;       
        end
    end
end
%