function dx = LorenzLorenzODENetdifR(x,a1,eps)
%INPUT
%-x: the initial conditions for the whole system in an array. 
%-a1: array of rho values for each Lorenz in the system. The rest of coefficents of the lorenz
%attractor are set to \simga=10 and \beta=8/3
%-eps: Matrix of couplings strengths.
%OUTPUT
%-dx :an array of all components of the Lorenz attractor for all the nodes of the system

scon=0;
nodes = length(eps);

for col=1:nodes
    deg(col)=sum(eps(:,col));
end

for con=0:3:3*nodes-1
scon=scon+1;
dx1(:,scon) = [10*(x(2+con)-x(1+con))+(-x(1+con)*deg(scon)+x(1:3:3*nodes)*eps(:,scon));a1(scon)*x(1+con)-x(2+con)-x(1+con)*x(3+con);x(1+con)*x(2+con)-(8/3)*x(3+con)];

end
dx = reshape(dx1,[],1);