clear all
n1=1.5075;%芯層折射率
n2=1.5;%包層折射率
ri=9.8*10^(-6);%最大的半徑
%r2=2.5*10^(-5);
dn=(n1-n2)/n1;
for ii=1:101
    r1(ii)=(ri/100)*ii;
end
n=n2.*ones(101,1);
for ii=1:60
    n(ii,1)=n1*sqrt((1-2*dn*(r1(ii)/r1(60))));
end
%[x,y]=meshgrid(1:101,1:101);mesh(n);
%for i=1:101 for j=1:101 if r1(i,j)<=r2 n(i,j)=n2; end end end

k0=2*pi/(1.55*10^(-6));
m=1;
zh=20*10^(-8);
n0=1.504;
for ii=1:101
    z1(ii)=ii*zh;
end
%ro
ah=0.268;%半徑漸變參數
for jj=1:101
    r(jj)=ri-ah*z1(jj);
end
for jj=1:101
rh(jj)=r(jj)/100;
end%注意
for jj=1:101
for ii=1:101
    b(ii,jj)=(k0^2)*(rh(jj)^2)*(n(ii,1)^2-(n0)^2)-2*(2-m^2/2);
end
end%系數b
for jj=1:101
for ii=1:101
        f(ii,jj)=(rh(jj)^2)*(k0^2)*(n(ii,1)^2-(n0)^2)-(rh(jj)^2)*(m^2)/r(jj)^2-2;
end%系數f
end
for jj=1:101
a(jj)=4*(k0^2)*(n0^2)*(rh(jj)^2)/(1+i*k0*n0*zh);
end
c=(1+exp(i*m*pi))*(2-m^2/2);
d=(1-k0*n0*i*zh)/(1+k0*n0*i*zh);
e=zeros(101,1);
g=zeros(101,1);
for ii=2:101
    e(ii,1)=1-1/(2*(ii-1));
end
for ii=2:101
      g(ii,1)=1+1/(2*(ii-1));  
end%系數整理結束
 %x=0;%光開始輸入
 %h1=0.03; 
 %for jj=1:101
  %  u(1,jj)=exp(-((x+jj)*h1)^2);
 %end%光結束
  r0=2.9815*10^(-6);
  for ii=1:101
     r1(ii)=(ri/100)*(ii-1);
  end
 for jj=1:101
    u(1,jj)=0.25*10^6*r1(jj)*exp(-3.33*10^5*(r1(jj)-r0))^2;
 end

for m1=1:101%大程序開始
A=zeros(101,101);
A(1,1)=a(m1)+b(1,m1);
A(1,2)=c;
for ii=2:101
    A(ii,ii-1)=e(ii,1);
end
for ii=2:100
    A(ii,ii)=f(ii,m1)+a(m1);
end
for ii=2:99
    A(ii,ii+1)=g(ii,1);
end%左面的系數矩陣A結束
B=zeros(101,101);
B(1,1)=a(m1)+b(1,m1)*d;
B(1,2)=c*d;
for ii=2:101
    B(ii,ii-1)=e(ii,1)*d;
end
for ii=2:100
    B(ii,ii)=f(ii,m1)*d+a(m1);
end
for ii=2:99
    B(ii,ii+1)=g(ii,1)*d;
end%右面的系數矩陣B結束
%d1=ones(101,101);
%y=ones(101,101);
d1(:,m1)=B*u.';%轉置
  d2=ones(101,101);
for ii=1:101
    d2(m1,ii)=d1(ii,m1);%變成一行一行的
end
    %下面用追趕法
  h2(m1,1)=c./(a(m1)+b(1,m1));
    for jj=2:100
    h2(m1,jj)=g(jj,1)./(f(1,m1)+a(m1)-e(jj,1)*h2(m1,jj-1));
    end
    y(m1,1)=d2(m1,1)/(a(m1)+b(1,m1));
    for jj=2:101
    y(m1,jj)=(d1(jj,m1)-e(jj,1).*y(m1,jj-1))./(f(1,m1)+a(m1)-e(jj,1).*h2(m1,jj-1));%problem
    end
    z(m1,101)=y(m1,101);
    for jj=100:-1:1
    z(m1,jj)=y(m1,jj)-h2(m1,jj).*z(m1,jj+1);
    end
    for jj=1:101
        u(1,jj)=z(m1,jj);
    end
    for jj=1:101
    q(m1,jj)=z(m1,jj);
    end
    end%最大循環的end
    [x,y]=meshgrid(0:1:100,0:zh:99*zh);
   % mesh(abs(q));
    mesh(abs(q));
    surf(abs(q));
    shading interp