%% 可以自動辨別前向跳變點，給出前向抽運閾值。n=7.88*10^20;m=8*10^19;segmaTm12=1*10^-23;segmaTm23=3.54*10^-21;wTm21=227;wTm31=2000;wTm32=304;wHo21=170;alafa=0;alafaTmHo=0;gamaTmHo=1*10^-17;gamaHoTm=3.9*10^-19;beta=3.1*10^-17;c=0.5;h=6.63*10^-34;mu=3.0*10^8/(1.06*10^-6);X0=[n;0;m];                 % difine the initial value for iterationi=0;Pn=40000;            % intensity range givendh=200;              % space step in iterationNum=Pn/dh+2;i=0;i_f=Num;                    % used to identify the divergent point                 N1f=zeros(Num-1,1);N2f=zeros(Num-1,1);M1f=zeros(Num-1,1);      % matrix for calculation results storage%% circling calculation by "for"for P=0:dh:Pn     % threshold maximum point        if (i+1)<i_f     % Jump out of the circuling calculation at divergent point        i=i+1;         % circulating number        % difine function group matrix and Jacobi matrix          yita=P/(h*mu);            f1=sym('yita*segmaTm23*n2-beta*n1*(n-n1-n2)-(wTm32+wTm31)*(n-n1-n2)+alafa*n2^2');        f2=sym('yita*segmaTm12*n1-yita*segmaTm23*n2+2*beta*n1*(n-n1-n2)+wTm32*(n-n1-n2)-wTm21*n2-gamaTmHo*n2*m1+gamaHoTm*n1*(m-m1)-2*alafa*n2^2-alafaTmHo*n2*(m-m1)+c*alafaTmHo*n2*(m-m1)');        f3=sym('-wHo21*(m-m1)+gamaTmHo*n2*m1-gamaHoTm*n1*(m-m1)-c*alafaTmHo*n2*(m-m1)');              FM=[f1;f2;f3];         DF11=diff(f1,'n1');DF12=diff(f1,'n2');DF13=diff(f1,'m1');        DF21=diff(f2,'n1');DF22=diff(f2,'n2');DF23=diff(f2,'m1');        DF31=diff(f3,'n1');DF32=diff(f3,'n2');DF33=diff(f3,'m1');            DFJ=[DF11,DF12,DF13;                                     DF21,DF22,DF23;             DF31,DF32,DF33];             n1=X0(1);        n2=X0(2);        m1=X0(3);     % initial value of levels concentration          F=subs(FM);        DF=subs(DFJ);        DFF=DF\F;      % with real matrix elements        X=X0-DFF;        a=norm(X-X0);        b=norm(X0);        c=a/b;        d=0;                  % Newton iteration for every specific "P"        while c>10^-2&d<1             % judge at the given accuracy with relative error and identify the convergence or divergence of iteration sequence                         X0=X;              X1=X0;             n1=X0(1);             n2=X0(2);             m1=X0(3);         %former value as the input value of later calculation             F=subs(FM);             DF=subs(DFJ);             DFF=DF\F;                       X=X0-DFF;             X2=X;             d1=norm(X2-X1);         % former space of the iteration sequence                           X0=X;              n1=X0(1);             n2=X0(2);             m1=X0(3);            %former value as the input value of later calculation             F=subs(FM);             DF=subs(DFJ);             DFF=DF\F;                       X=X0-DFF;             X3=X;             d2=norm(X3-X2);       % later space of the iteration sequence                   d=norm(d2/d1);        % the ratio of later and former sequence             a=norm(X-X0);             b=norm(X0);             c=a/b;           end                if c<=10^-2           N1f(i)=X(1);           N2f(i)=X(2);           M1f(i)=X(3);     % producing matrix elements for every levle concentration           X0=X;           % former value as the input value of later calculation       else           disp('>>****************************************************************')           disp('     A divergent sequence occurs as increasing pumping density ! ')           disp('     The below corresponds to the first divergent point ')                             disp('>>****************************************************************')           i_f=i;           Pmax=P-dh;           IfP=[i_f Pmax]              % show the first divergent point and forward pumping threshold       end          endendN3f=n*ones(Num-1,1)-N1f-N2f;M2f=m*ones(Num-1,1)-M1f;N1f;N2f;N3f;M1f;M2f; %% 可以自動辨別逆向跳變點，給出逆向抽運閾值。 n1P=1.78*10^20;n2P=6*10^20;m1P=0.01*10^19;              % given initial values of level concentration for power intensity much more than thresholdX0=[n1P;n2P;m1P];             % difine the initial value for iterationPnb=40000;i=Pnb/dh+2;i_b=0;N1b=zeros(Num-1,1);N2b=zeros(Num-1,1);M1b=zeros(Num-1,1);      % matrix for calculation results storage%% circling calculation by "for"for P=Pnb:-dh:0        if (i-1)>=i_b               i=i-1;       % difine function group matrix and Jacobi matrix         yita=P/(h*mu);           f1=sym('yita*segmaTm23*n2-beta*n1*(n-n1-n2)-(wTm32+wTm31)*(n-n1-n2)+alafa*n2^2');       f2=sym('yita*segmaTm12*n1-yita*segmaTm23*n2+2*beta*n1*(n-n1-n2)+wTm32*(n-n1-n2)-wTm21*n2-gamaTmHo*n2*m1+gamaHoTm*n1*(m-m1)-2*alafa*n2^2-alafaTmHo*n2*(m-m1)+c*alafaTmHo*n2*(m-m1)');       f3=sym('-wHo21*(m-m1)+gamaTmHo*n2*m1-gamaHoTm*n1*(m-m1)-c*alafaTmHo*n2*(m-m1)');             FM=[f1;f2;f3];        DF11=diff(f1,'n1');DF12=diff(f1,'n2');DF13=diff(f1,'m1');       DF21=diff(f2,'n1');DF22=diff(f2,'n2');DF23=diff(f2,'m1');       DF31=diff(f3,'n1');DF32=diff(f3,'n2');DF33=diff(f3,'m1');           DFJ=[DF11,DF12,DF13;                                    DF21,DF22,DF23;            DF31,DF32,DF33];              n1=X0(1);       n2=X0(2);       m1=X0(3);    % initial value of levels concentration         F=subs(FM);       DF=subs(DFJ);       DFF=DF\F; % with real matrix elements       X=X0-DFF;       a=norm(X-X0);       b=norm(X0);       c=a/b;       d=0;         % Newton iteration for every specific "P"      while c>10^(-2)&d<1             % judge at the given accuracy with relative error            X0=X;            X1=X0;            n1=X0(1);            n2=X0(2);            m1=X0(3); %former value as the input value of later calculation            F=subs(FM);            DF=subs(DFJ);            DFF=DF\F;                      X=X0-DFF;            X2=X;            d1=norm(X2-X1);         % former space of the iteration sequence                            X0=X;             n1=X0(1);            n2=X0(2);            m1=X0(3);            %former value as the input value of later calculation            F=subs(FM);            DF=subs(DFJ);            DFF=DF\F;                      X=X0-DFF;            X3=X;            d2=norm(X3-X2);       % later space of the iteration sequence                  d=norm(d2/d1);        % the ratio of later and former sequence            a=norm(X-X0);            b=norm(X0);            c=a/b;           end                if c<=10^-2            N1b(i)=X(1);            N2b(i)=X(2);            M1b(i)=X(3);     % producing matrix elements for every levle concentration            X0=X;           % former value as the input value of later calculation         else            disp('>>****************************************************************')            disp('     A divergent sequence occurs as decreasing pumping density ! ')            disp('     The below corresponds to the first divergent point ')                              disp('>>****************************************************************')            i_b=i;            Pmin=P+dh;            IbP=[i_b Pmin]              % show the first divergent point and forward pumping threshold         end    endendN3b=n*ones(Num-1,1)-N1b-N2b;M2b=m*ones(Num-1,1)-M1b;N1b;N2b;N3b;M1b;M2b; % 自動整合前向和逆向數據點a=1:i_b;N1b(a)=N1f(a);N2b(a)=N2f(a);M1b(a)=M1f(a);M2b(a)=M2f(a);b=i_f:Num-1;N1f(b)=N1b(b);N2f(b)=N2b(b);M1f(b)=M1b(b);M2f(b)=M2b(b);% show figuresP=0:dh:Pn;figure;hold on;plot(P,N1f,'k');plot(P,N2f,'r');plot(P,10*M1f,'b');plot(P,10*M2f,'g');plot(P,N1b,'k');plot(P,N2b,'r');plot(P,10*M1b,'b');plot(P,10*M2b,'g'); % 保存數據save c:\matlab6p1\work\SSL-IOB\fig3data-f N1f N2f N3f M1f M2f   % save the forward results calculatedsave c:\matlab6p1\work\SSL-IOB\fig3data-b N1b N2b N3b M1b M2b   % save the backward results calculated