计算电磁学王秉中矩形波导基膜和第一高次模色散特性曲线

clear;clc;close all;
sz_X=41;sz_Y=23;
sz_M=sz_X-2;sz_N=sz_Y-2;
h=0.254;
% sz_X=81;sz_Y=45;
% sz_M=sz_X-2;sz_N=sz_Y-2;
% h=0.127
% sz_X=21;sz_Y=12;
% sz_M=sz_X-2;sz_N=sz_Y-2;
% h=0.508;
%TM波计算
DM=diag（4*ones（sz_M1）0）+diag（-ones（sz_M-11）1）+diag（-ones（sz_M-11）-1）;
for i=1:sz_N
    D1（::i）=DM;
end
 sz_D1 = size（D1）;
c = mat2cell（reshape（D1sz_D1（1）[]）sz_D1（1）repmat（sz_D1（2）1sz_D1（3）））;
 K1=feval（@blkdiagc{:}）;
 KM=K1+diag（-ones（1sz_M*sz_N-sz_M）sz_M）+diag（-ones（1sz_M*sz_N-sz_M）-sz_M）;
 [v_TMd_TM]=eig（KM）;
 [dx_TMdy_TM]=find（d_TM>0）;
 l_dxTM=length（dx_TM）;
 for i=1:l_dxTM
     d2_TM（i）=d_TM（dx_TM（i）dy_TM（i））;
 end
 for i=1:l_dxTM-1
     minIndex=i;
     for j=i+1:l_dxTM
         if（d2_TM（j）             minIndex=j;
         end
     end
     temp=d2_TM（i）;
     d2_TM（i）=d2_TM（minIndex）;
     d2_TM（minIndex）=temp;
 end
 kc_TM=sqrt（d2_TM）/h;
 %TE波计算
 DE=diag（4*ones（sz_X1）0）+diag（-ones（sz_X-11）1）+diag（-ones（sz_X-11）-1）;
 DE（12）=-2;DE（sz_Xsz_X-1）=-2;
for i=1:sz_Y
    D2（::i）=DE;
end
 sz_D2 = size（D2）;
c = mat2cell（reshape（D2sz_D2（1）[]）sz_D2（1）repmat（sz_D2（2）1sz_D2（3）））;
 K1=feval（@blkdiagc{:}）;
 KE=K1+diag（-1*ones（1sz_X*sz_Y-sz_X）sz_X）+diag（-1*ones（1sz_X*sz_Y-sz_X）-sz_X）;
 for i=1:sz_X
     KE（ii+sz_X）=-2;
     KE（i+sz_X*（sz_Y-1）i+sz_X*（sz_Y-2））=-2; %KE矩阵
 end
 [v_TEd_TE]=eig（KE）;
 [dx_TEdy_TE]=find（d_TE>1e-6）;
 l_dxTE=length（dx_TE）;
 for i=1:l_dxTE
     d2_TE（i）=d_TE（dx_TE（i）dy_TE（i））;
 end
 for i=1:l_dxTE-1
     minIndex=i;
     for j=i+1:l_dxTE
         if（d2_TE（j）