function dxdt=mlcnRHS(t,x,p)

  dxdt=zeros(3*p.n,1);
  for i=0:(p.n-1)  
    dxdt(3*i+1)=x(3*i+2)/p.L-p.r0*x(3*i+1)/p.L;
    dxdt(3*i+2)=(-x(3*i+1)-p.alpha*x(3*i+2)*x(3*i+3))/p.C;
    if p.n>1
      j=mod(i-1,p.n);
      dxdt(3*i+2)=dxdt(3*i+2)+p.G*(x(3*j+2)-x(3*i+2))/p.C;
    end
    if p.n>2
      j=mod(i+1,p.n);
      dxdt(3*i+2)=dxdt(3*i+2)+p.G*(x(3*j+2)-x(3*i+2))/p.C;
    end
    dxdt(3*i+3)=x(3*i+2)*x(3*i+2)-p.beta;
  end
