Help for new type of ode45 question: coupled equations but each with large summations on dependent variables

Hello everyone,

I have seen all ode45 answers but no case seems to be like the one below in which each equation has summations comprising of other 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.

    clc
    clear all
%     dz= 400; nz=200;
%     L= nz*dz;
%     zspan = linspace(0,L,nz+1)  %was trying to some other way of zspan

    ws=[1520:10:1560]*1e-9 ;    F=3e8./ws;  %wavelength and frequency of the channels
    channels=numel(F);   num_span=1;
    P_dbm= -3.5* ones(channels, num_span);  %input power of each channel
  
    p0= 10.^(P_dbm./10)*0.001;  %convert input signal power to Watts
    pp1=300e-3 ;        pp2=310e-3 ;   %300 mW and 310 mW initial pump powers
    p_in=[p0;pp1;pp2];   %initial condition vector having signal powers for 5 channels+ two pump powers
    zspan = [0:40]*1000   %for 40 kms   THESE BIG NOS. ARE CAUSING PROBLEMS
    
    
    [z,p]= ode45(@lorenz,zspan,p_in);  %the ode calling
    p_in_dBm=(10*log10(p))+30;   %back to dBm
      
    %%
    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