How to find mode shape

this is my code for one shaft 0 disk. but I don't know how to plot mode shape I can't find b.

 nel1=14;                    
 nnel=2;                    % 
ndof=4;                    % 
phase1=0;                  % 
nb1=2;                     % 
snode1=nel1+1;             % 
nnode1=(nel1+1);           %
bpo1=[0,1.4];              % 
edof=nnel*ndof;            % 
sdof=snode1*ndof;          % 
L1=1.4;                %
Es1=2e11;         %　
v1=0.3;             %　
Ros1=7800;           %　
Gs1=Es1/(2*(1+v1));  %　
G1=Gs1;                               
rr1=zeros(2,nel1);   %　
rr1(1,:)=0.05;      %　
rr1(2,:)=0;          %　

 MM1=zeros(sdof,sdof);        % initialization of M system matrix
CC1=zeros(sdof,sdof);        % initialization of C system matrix
KK1=zeros(sdof,sdof);        % initialization of K system matrix
GG1=zeros(sdof,sdof)
 gcoord1=zeros(nnode1,1);
 for i=1:nnode1;
      gcoord1(i)=gcoord1(i)+(i-1)/nel1*L1;
 end

 for i=1:nel1;
    for j=1:nnel;
    nodes1(i,j)=i+j-1;
    end
end
 for iel=1:nel1;
    nl1=nodes1(iel,1); nr1=nodes1(iel,2); 
    L21=gcoord1(nr1)-gcoord1(nl1);        
    L11=gcoord1(nl1);                  
    rso1=rr1(1,iel);
    rsi1=rr1(2,iel);
    m1=rsi1/rso1;                 
    kpran1=(6*(1+v1)*(1+m1^2)^2)/((7+6*v1)*(1+m1^2)^2+(20+12*v1)*m1^2);  
    o1=(12*Es1*(1/4*pi*(rso1^4-rsi1^4)))/(kpran1*Gs1*(pi*(rso1^2-rsi1^2))*L21^2);
    As1=pi*(rso1^2-rsi1^2);
    Isd1=1/4*pi*(rso1^4-rsi1^4);
    Isp1=1/2*pi*(rso1^4-rsi1^4);
    [m1,k1,g1]=residualMCK1(o1,L1,Ros1,L11,L21,As1,Isp1,Isd1,Es1,kpran1,Gs1); 
    mm1=m1;
    kk1=k1;
    gg1=g1;

 for i=1:edof 
 start1 = (iel-1)*(nnel-1)*ndof;
      ii=start1+i;
      for j=1:edof
         jj=start1+j;
         MM1(ii,jj)=MM1(ii,jj)+double(mm1(i,j));
         KK1(ii,jj)=KK1(ii,jj)+double(kk1(i,j));     
         GG1(ii,jj)=GG1(ii,jj)+double(gg1(i,j)); 
     end
  end
end
 KKB1=zeros(sdof,sdof);
 KKB2=zeros(sdof,sdof);
 CCB1=zeros(sdof,sdof);
 CCB2=zeros(sdof,sdof);

 Cyy1=0;
 Cyy2=0;
 Cyz1=0;
 Czy1=0;
 Czz1=0;
 Czz2=0;
 Kyy1=1e12;
 Kyy2=1e12;
 Kyz1=0;
 Kzy1=0;
 Kzz1=1e12;
 Kzz2=1e12;

 Kb1=[Kyy1 Kyz1;Kzy1 Kzz1];
 Cb1=[Cyy1 Cyz1;Czy1 Czz1];
 Kb2=[Kyy2 Kyz1;Kzy1 Kzz2];
 Cb2=[Cyy2 Cyz1;Czy1 Czz2];
 for ib1=1:nb1;
     start1=bpo1(ib1)/L1*nel1*ndof;
    for i=1:2;
        ii=start1+i;
        for j=1:2;
            jj=start1+j;
            KKB1(ii,jj)=KKB1(ii,jj)+Kb1(i,j);     %  combine bearing to global
            CCB1(ii,jj)=CCB1(ii,jj)+Cb1(i,j);
            KKB2(ii,jj)=KKB2(ii,jj)+Kb2(i,j);     %  combine bearing to global
            CCB2(ii,jj)=CCB2(ii,jj)+Cb2(i,j);
       end
    end
end
 MMMM=MM1;
 CCCC1=CC1+CCB1;
 CCCC2=CC1+CCB2;
 KKKK1=KK1+KKB1;
 KKKK2=KK1+KKB2;
 GGGG=GG1;
 o=zeros(sdof,sdof);
 n=0;
 XXYY=zeros(1,1);

    Oma1=0;
    n=n+1;
    XXYY(n)=Oma1;%*60/2/pi;
    E=[o,-MMMM;MMMM,(CCCC1)+(Oma1*GGGG)];
    F=[MMMM,o;o,KKKK1];
    EF=-inv(E)*F;
    w=sort(abs(imag(eig(EF))));
    wwwa=abs(imag(eig(EF)));
    [C,D]=eig(EF);
    K=imag(C);
    R=real(C);
    ww1(:,n)=w;%*60/2/pi;

 A0=imag(D);
 for k=1:sdof*2
     A1(k,1)=A0(k,k);
 end

 for i=1:2:sdof*2-1
  QQQ(:,i)=C(:,i)*D(i,i)+C(:,i+1)*D(i+1,i+1); 

 end
 QQQ(:,sdof*2)=C(:,sdof*2)*D(sdof*2,sdof*2);

 QQQV1=QQQ(sdof+1:ndof:sdof*2,:);
 QQQW1=QQQ(sdof+2:ndof:sdof*2,:);
 QQQB1=QQQ(sdof+3:ndof:sdof*2,:);
 QQQR1=QQQ(sdof+4:ndof:sdof*2,:);

 QQQV2=QQQ(sdof*1.5+1:ndof:sdof*2,:);
 QQQW2=QQQ(sdof*1.5+2:ndof:sdof*2,:);
 QQQU2=QQQ(sdof*1.5+3:ndof:sdof*2,:);
 QQQB2=QQQ(sdof*1.5+4:ndof:sdof*2,:);
 QQQR2=QQQ(sdof*1.5+5:ndof:sdof*2,:);
 QQQA2=QQQ(sdof*1.5+6:ndof:sdof*2,:);

% plot mode shape % figure(1) %coupled b=117; plot(1:snode1,QQQV1(:,b),1:snode1,QQQW1(:,b),'linewidth',0.5); title('Mode 1 = 633.65 (rad/s)') legend('Driving-V','Driving-W') ;

 figure(5)
 %coupled
 b=31;
 plot(1:snode1,QQQV1(:,b),1:snode1,QQQW1(:,b),'linewidth',0.5);
 title('Mode 2 = 2550 (rad/s)')
 legend('Driving-V','Driving-W') ;

[KSSV]

We cannot run your code...few functions are missing.