New type of ode45 question: coupled equations but with large summations on dependent variables

Question

Jaya 2021-4-25

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/812765-new-type-of-ode45-question-coupled-equations-but-with-large-summations-on-dependent-variables

I am trying to replicate the following set of equations. i am getting wrong values. Can anyone see the code and please tell me if the method of ode45 implementation is correct? If correct, then it would mean there is something wrong with the input values that is causing wrong results for me. I have commented everything for better understanding.

% Here, I define the zspan and p_initial and call the following function
    
    function dpdz= lorenz(z,p)
    %%
    ws=[1520:10:1560]*1e-9;    F=3e8./ws;
    f_p1= 3e8/1427e-9;      f_p2=3e8/1495e-9;
    gr=.4/1e3;     %Raman gain spectrum
    alphaa=[0.1832    0.1932    0.1142    0.0746    0.1781    0.1432    0.1969]./4.343/1e3 ;
    %attenuation coefficient encountered by each frequency channel
    channels=numel(F);
    m=length(p);
    for i=1: channels      %For each Pi. The first equation
        A=0; B=0;   f_i=F(i);
        for k= i+1:channels     %I separately wrote the components of each eqn. and then added them accordingly
            A= A-((f_i/F(k))*p(k)*p(i))*gr;
        end
        for k=1:(i-1)
            B= B+(p(k)*p(i)*gr);
        end
        alpha_i=alphaa(i);
        C= p(i)*alpha_i;
        D= p(m-1)*p(i)*gr;
        E=p(m)*p(i)*gr;
        dpdz(i)= A+B-C+D+E;    %this gives one equation. As 5 channels, dpdz(1) to dpdz(5), five eqns.
    end   %end of the every channel P loop
    %=======================================================%
    %for the pump 1 eqn.
    A=0; B=0; C=0;
    for k=1:channels
        A= A+(((f_p1/F(k))*p(k)*p(m-1))*gr) ;
    end
    B= -p(m)*p(m-1)*gr;     
    C=p(m-1)*alpha_i; 
    dpdz(m-1)= A+B+C;
    %=======================================================%
    %for the pump 2 eqn.
    A=0; B=0; C=0;
    for k=1:channels
        A= A+ (((f_p2/F(k))*p(k)*p(m))*gr) ;
    end
    B= (f_p2/f_p1)*p(m)*p(m-1)*gr;   %need to add the gr term later
    C=p(m)*alpha_i;  %later you will have to change the alpha notation
    dpdz(m)= A+B+C;
    dpdz=dpdz';
    %=======================================================%
      end