EC=1^2;
ES=3.8^2;
ER=5.20^2;
landa=1.55;
k0=2*pi/landa;
hx=0.05;
hy=0.05;
hx=0.05;
xlo=20;
xhi=35;
ylo=15;
yhi=45;
NX=50;
NY=59;
prm=EC*ones(NX,NY);
prm([1:xlo-1],

=ES*ones(xlo-1,NY);
prm([xlo+1:xhi-1],[ylo+1:yhi-1])=ER*ones(xhi-xlo-1,yhi-ylo-1);
prm(xlo)=0.5*(prm(xlo-1)+prm(xlo+1));
prm(xhi,[ylo+1:yhi-1])=0.5*(ER+EC)*ones(1,yhi-ylo-1);
prm([xlo+1:xhi-1],[ylo,yhi])=0.5*(ER+EC)*ones(xhi-xlo-1,2);
prm(xlo,[ylo,yhi])=0.25*(2*ES+EC+ER);
prm(xhi,[ylo,yhi])=0.25*(3*EC+ER);
functionH=helmholtz(hx,hy,d)
[NX,NY]=size(d);
N=NX*NY;
ihx2=1/hx/hx;
ihy2=1/hy/hy;
md=-2*(ihx2+ihy2)*ones(N,1)+reshape(d,NX*NY,1);
xd=ihx2*ones(N,1);
yd=ihy2*ones(N,1);
H=spdiags([yd,xd,md,xd,yd],[-NX,-1,0,1,NX],N,N);
for n=NX:NX:N-NX
H(n,n+1)=0;
H(n+1,n)=0;
end