clear
a=3.6;
b=4;
d=5.8;
k=(4*pi)/3;
w=12.5716e08;
e0=8.85e-12;
y=1.35;
x=[linspace(10,50,5)];
z=[linspace(10,50,5)];
x1=0.45;
y1=1.4;
z1=1.37;
e0n=2;
m = [linspace(4,100,97)];
n = [linspace(4,100,97)];
kx=(m*pi)/a;
ky=(n*pi)/b;
beta=sqrt(kx.^2+ky.^2-k^2);
E0=1/(w*e0);
E1=(2*e0n*sin(kx*x1).*cos(ky*y1))/(a*b*beta.*sin(beta*d));
E2=sum(E1);
E4=zeros(6);
E=zeros(6);
for x=10:10:50
for z=10:10:50
E3=(kx.*ky.*cos(kx.*x).*sin(ky*y).*sin(beta*z1).*sin(beta.*(d-z)))+((ky.^2-k^2).*sin(kx.*x).*cos(ky*y).*sin(beta*z1).*sin(beta.*(d-z)))+(ky.*beta.*sin(kx.*x).*sin(ky*y).*(-sin(beta*z1)).*cos(beta.*(d-z)));
E4(x,z)=sum(E3);
end
end
E(x,z)=E0*E2.*E4(x,z);
plot3(z,x,E)

(please let me know the corrections to make.thanks)