how to add a loss term to my function?

my code currently looks like:

 function Vdot = FEL_equationFifthYear(z,V);
global N p0 a0 phi0 Re ku K gamma0 rho;

x = cos(V(1:N,1) + V(2*N+2,1));
y = sin(V(1:N,1) + V(2*N+2,1));

Vdot = [1:2*N+1,1]';
Vdot(1:N,1) = V((N+1):2*N,1);
Vdot(N+1:2*N,1) = ((-2*V(2*N+1,1)*cos(V(1:N,1) + V(2*N+2,1))))-((2/3)*Re*(((1/N)*sum(gamma0))*(1+rho*((1/N)*sum(V((N+1):2*N,1)))))^2*ku^2*K^2);;
Vdot(2*N+1,1) = (1/N)*sum(x);
Vdot(2*N+2,1) = -(1/V(2*N+1,1))*((1/N)*sum(y));
end

 clear all;
global N p0 a0 phi0 rho gamma0 Re ku K loss %N = number of electrons, 

 N=50; %Number of electrons
p0=0; %Pj=P0=constant. Used 0 for ease of use but from worked equations Pj=5.67e7, is this 0 up to 5.67e7. working in jotter
a0=10^-3;%A(z=0)=A0<<1, so 0.0001 was chosen 
phi0 = 0;
rho = 0;
gamma0 = 0;
Re = 2.818e-15;
ku = 1;
K = 3.71;
loss = -(((2/3)*Re*(((1/N)*sum(gamma0))*(1+rho*((1/N)*sum(p0))))^2*ku^2*K^2));

 V0 = [1:2*N+2]'; %Vector V from 1 to 2*N+2, apostrophe used to invert the matrice. This shows the complete vector V
V0(1:N,1) = linspace(0,((N-1)/N)*2*pi,N); %vector V from 1:N evenly spaced out between points X1 and X2 with N number of points
V0((N+1):2*N,1) = p0 + loss; %Vector V from point N+1 to 2N is equal to p0
V0(2*N+1,1) = a0;
V0(2*N+2,1) = phi0;

 zspan = [0 200];

 z0=[0:101]*0;
options=odeset();
[z,V] = ode45('FEL_equationFifthYear',zspan,V0,options); 

 %V(1:N,1)= pj;
%V(N+1:2*N,1) = -(A*exp(ithetaj)+c.c.)+loss;
%V(2N+1,1) = <exp(itheta)>
%plot(z,(1/N)*sum(V(N+1:2*N,1)));
ylim([-10 10]);

 semilogy(z, V(:,2*N+1).^2);
%title('graph of Ln(A)^2 Vs z(bar)')
%xlabel('z(bar)')
ylabel('Ln(A)^2')

The problem comes when I want to change my loss term: loss = -(((2/3)*Re*(((1/N)*sum(gamma0))*(1+rho*((1/N)*sum(p0))))^2*ku^2*K^2)); from p0 to pj. This pj term is introduced after the ode function as V(1:N,1). how would I go about adding this term to my code?