function [zeta,s,U]=diffraction(m,n)
dzeta=1/m;
ds=6/n;
ta=dzeta/(ds*ds);
zeta=zeros(m+1,1);
s=zeros(n-1,1);
U=zeros(m+1,n-1);
b=zeros(n-1,1);
ll=zeros(n-1,1);
ul=zeros(n-1,1);
zeta=[0:dzeta:1];
s=[-(n-2)/2*ds:ds
n-2)/2*ds]';
%r=3;
%thet0=4.5*10^-3;
f=inline('sqrt(0.1/sqrt(pi)/(4*10^-3)*exp(-(2*10^-3).^2/(4*10^-3).^2))*exp(-s.^2/(0.69*0.69))','s');
for j=1:n-1
U(1,j)=f(s(j));
end;
ne=2.3;
lambda0=0.5*10^-6;
ke=2*pi/lambda0*ne;
thet=-2*10^-3;
x0=18.265*10^-6;
alpha=thet*ke*x0;
g0=alpha/2*ds*ta;
g1=1/4*sqrt(-1)*ta;
for k=1:n-2
A(k,k+1)=g0-g1;
A(k+1,k)=-g1;
A(k,k)=-g0+1+2*g1;
end
A(1,1)=A(1,1)-g1;
A(n-1,n-1)=1+g1;
% 分解矩阵A
ul=A(1,1);
for l=2:n-1;
ll(l)=A(l,l-1)/ul(l-1);
ul(l)=A(l,l)-A(l-1,l)*ll(l);
end
for k=1:m;
b(1)=(1+g0-g1)*U(k,1)+(-g0+g1)*U(k,2);
for j=2:n-2;
b(j)=g1*U(k,j-1)+(1+g0-2*g1)*U(k,j)+(-g0+g1)*U(k,j+1);
end
b(n-1)=g1*U(k,n-2)+(1-g1)*U(k,n-1);
%用追赶法求解线性方程组AU=b
y(1)=b(1);
for l=2:n-1
y(l)=b(l)-ll(l)*y(l-1);
end;
U(k+1,n-1)=y(n-1)/ul(n-1);
for l=n-2
-1):1
U(k+1,l)=(y(l)+(-g0+g1)*U(k+1,l+1))/ul(l);
end;
end
mesh(s,zeta,abs(U).^2)
xlabel('s')
ylabel('\xi')
zlabel('|U|.^2')