clear
format long

global delta Neff delta_neff_average tlamda Pg Pe epsrate

epsrate=1.0e-10;

Pe=0.220;

L=4e+6; % 4mm

Neff=1.45;
lamda_b=1550;
Pg=lamda_b/(2*Neff);
delta_neff_average=2e-4;
lamda=(1548:0.01:1552);
Mlamda=length(lamda);
step=lamda(2)-lamda(1);

zspan=[L,0];
rz0=0;
opts=odeset('reltol',1.0e-6);

for ilam=1:Mlamda
    
fprintf(1,'ilam= %d \n',ilam); 
 
lamda(ilam)=lamda(1)+(ilam-1)*step;

tlamda=lamda(ilam);

delta=2.0*pi*Neff*(1.0/lamda(ilam)-1.0/lamda_b);
  
[z,r]=ode45('drdz',zspan,rz0,opts);

rr(ilam)=r(end);

refl(ilam,1)=lamda(ilam);
refl(ilam,2)=abs(rr(ilam))^2;

end