%計算光纖激光器的正反向功率與位置的關系
clear
z0=0;

global input
   input=3;
for i=1:20
l(i)=(i)*50
z1=l(i);
  a01=[2,2/0.99];
  j=1;t=1;
  while j==1 
     t=t+1;
      a02=a01+[1e-5,1e-5/0.99];
    [z,p]=ode45('laser18',[z0 z1],a01); 
    n=length(z);
    p1=p(:,1);                     %正向功率
    p2=p(:,2);                     %反向功率
  
   h1=p2(n)-0.04*p1(n);
   [x,w]=ode45('laser18',[z0 z1],a02);
   m=length(x);
   w1=w(:,1);
   w2=w(:,2);
 
   h2=w2(m)-0.04*w1(m);
    if abs(h1)<1e-5  
        break;
    else 
       a01(1)=a01(1)-h1/((h2-h1)/1e-5);
       a01=[a01(1) a01(1)/0.99];
   end
end
out(i)=p1(n)-p2(n);
end

hold on
plot(l/100,out,'-')
%plot(z/100,p1,'-')
hold on
 
xlabel('Fiber length/m')
%plot(r1,pw1,'>')
ylabel('Output power /W')  
%legend('theory','experiment')
%p1(n)-p2(n)