matlab计算闪耀光栅TE波衍射效率

利用严格耦合波理论计算闪耀光栅TE波，TM波的衍射效率，并用matlab建模计算。下面是我程序不知道错哪了

你好，请问可以把报出错误的代码段截个图吗。然后就是如果没有报错的话，可不可以讲一讲你的预期和现在的不同之处，最好提供一下源代码，粘贴如以下格式：
（点右上方</>两行三点之间插入你的代码）

clear
n1 = 1.5;

clear
n1=1.52;%入射介质折射率
ep1=n1^2;
n2=1.52;%光栅折射率
ep2=n2^2;
n3=1;
ep3=n3^2;
lamd=5501e-9;%入射光波长
Lamd=5411e-9;%光栅周期
d=2081e-9;
N=20;%层数
y=zeros(1,N);%中位点
fg=linspace(0,d,N+1);
for m=1:1:N
y(m)=(fg(m)+fg(m+1))/2;
end
x=Lamd-yLamd/d;%x坐标
f=x/Lamd;%占空比
dn=d/N;%层厚
M=3;%衍射级数目
i=(-(M-1)/2:(M-1)/2)';%衍射级次
p=(M+1)/2;%0级衍射位置
theta=0;%入射角
k0=2pi/lamd;
K=2pi/Lamd;
kxi=k0n1sind(theta)+iK;
k1zi=conj(k0sqrt(n1^2-(kxi/k0).^2));
k3zi=conj(k0sqrt(n3^2-(kxi/k0).^2));
Kx=diag(kxi/k0);
Y1=diag(k1zi/k0);
Y2=diag(k3zi/k0);
h=(1-M:M-1)';
El=zeros(M,M,N);
Al=zeros(M,M,N);
Q=zeros(M,M,N);
W=zeros(M,M,N);
Q1=zeros(M,M,N);
V=zeros(M,M,N);
q=zeros(M,1,N);
X=zeros(M,M,N);
B=zeros(2M,2M,N);
D=zeros(2M,2M,N);
F=zeros(2M,2M,N);
for l=N:-1:1
ep=(n2^2-n3^2).sin(pif.h)./(pi.h);
ep(M,:)=n2^2f+n3^2(1-f);
El(:,:,l)=toeplitz(flipud(ep(1:M,l)),ep(M:2M-1,l));
Al(:,:,l)=Kx^2-El(:,:,l);
[W(:,:,l),Q(:,:,l)]=eig(Al(:,:,l),'nobalance');
Q1(:,:,l)=sqrt(Q(:,:,l));
q(:,1,l)=diag(Q1(:,:,l));
V(:,:,l)=W(:,:,l)Q1(:,:,l);
X(:,:,l)=diag(exp(k0q(:,1,l)*dn));
B(:,:,l)=[W(:,:,l),W(:,:,l)*X(:,:,l);-V(:,:,l), V(:,:,l)*X(:,:,l)];
D(:,:,l)=[W(:,:,l)*X(:,:,l),W(:,:,l);-V(:,:,l)*X(:,:,l),V(:,:,l)];
F(:,:,l)=B(:,:,l)/(D(:,:,l));

end
G=eye(2M);
for j=1:N
G=GF(:,:,j);
end
I=eye(M);
G1=G*[I;1iY2];
G2=[I;-1iY1];
C1=1iY1G1(1:M,:);
C2=G1(M+1:2M,:);
C=C1+C2;
L=[zeros(p-1,1);1i(k1zi(p)/k0+n1cosd(theta));zeros(p-1,1)];
T=C\L;
L2=[zeros(p-1,1);1;zeros(p-1,1);zeros(p-1,1);1in1cosd(theta);zeros(p-1,1)];
R=G2(G1T-L2);
DEri=R.conj(R).real(k1zi/(k0n1cosd(theta)));
DEti=T.conj(T). real(k3zi/(k0n1cosd(theta)));

下面是推导公式

clear
n1=1.52;%入射介质折射率
ep1=n1^2;
n2=1.52;%光栅折射率
ep2=n2^2;
n3=1;
ep3=n3^2;
lamd=550*1e-9;%入射光波长
Lamd=541*1e-9;%光栅周期
d=208*1e-9;
N=20;%层数
y=zeros(1,N);%中位点
fg=linspace(0,d,N+1);
for m=1:1:N 
    y(m)=(fg(m)+fg(m+1))/2;
end
x=Lamd-y*Lamd/d;%x坐标
f=x/Lamd;%占空比
dn=d/N;%层厚
M=3;%衍射级数目
i=(-(M-1)/2:(M-1)/2)';%衍射级次
p=(M+1)/2;%0级衍射位置
theta=0;%入射角
k0=2*pi/lamd;
K=2*pi/Lamd;
kxi=k0*n1*sind(theta)+i*K;
k1zi=conj(k0*sqrt(n1^2-(kxi/k0).^2));
k3zi=conj(k0*sqrt(n3^2-(kxi/k0).^2));
Kx=diag(kxi/k0);
Y1=diag(k1zi/k0);
Y2=diag(k3zi/k0);
h=(1-M:M-1)';
El=zeros(M,M,N);
Al=zeros(M,M,N);
Q=zeros(M,M,N);
W=zeros(M,M,N);
Q1=zeros(M,M,N);
V=zeros(M,M,N);
q=zeros(M,1,N);
X=zeros(M,M,N);
B=zeros(2*M,2*M,N);
D=zeros(2*M,2*M,N);
F=zeros(2*M,2*M,N);
for l=N:-1:1
    ep=(n2^2-n3^2).*sin(pi*f.*h)./(pi.*h);
    ep(M,:)=n2^2*f+n3^2*(1-f);
    El(:,:,l)=toeplitz(flipud(ep(1:M,l)),ep(M:2*M-1,l));
    Al(:,:,l)=Kx^2-El(:,:,l);
    [W(:,:,l),Q(:,:,l)]=eig(Al(:,:,l),'nobalance');
    Q1(:,:,l)=sqrt(Q(:,:,l));
    q(:,1,l)=diag(Q1(:,:,l));
    V(:,:,l)=W(:,:,l)*Q1(:,:,l);
    X(:,:,l)=diag(exp(k0*q(:,1,l)*dn));
    B(:,:,l)=[W(:,:,l),W(:,:,l)*X(:,:,l);-V(:,:,l), V(:,:,l)*X(:,:,l)];
    D(:,:,l)=[W(:,:,l)*X(:,:,l),W(:,:,l);-V(:,:,l)*X(:,:,l),V(:,:,l)];
    F(:,:,l)=B(:,:,l)/(D(:,:,l));
    
end
G=eye(2*M);
for j=1:N
   G=G*F(:,:,j);
end
I=eye(M);
G1=G*[I;1i*Y2];
G2=[I;-1i*Y1];
C1=1i*Y1*G1(1:M,:);
C2=G1(M+1:2*M,:);
C=C1+C2;
L=[zeros(p-1,1);1i*(k1zi(p)/k0+n1*cosd(theta));zeros(p-1,1)];
T=C\L;
L2=[zeros(p-1,1);1;zeros(p-1,1);zeros(p-1,1);1i*n1*cosd(theta);zeros(p-1,1)];
R=G2\(G1*T-L2);
DEri=R.*conj(R).*real(k1zi/(k0*n1*cosd(theta)));
DEti=T.*conj(T).* real(k3zi/(k0*n1*cosd(theta)));