function [storex,store,phi] =  RLS_net2(N,z0,g,q,omega,phi,eta,nt,imax,imin,P,m,dt,sup,train,inp,win)
z = z0;

%% run the training procedure. 
for j = 1:nt 
r = (z.*(z>0)).^m;
x = phi'*r;
z = z + dt*(-z + omega*(g*r)+q*eta*x+win*inp(j,1));    
store(j,:) = r(1:10); 
storex(j,:) = x;




if train==1 
if mod(j,3)==1
  if j>imin
      if j < imax
        
        e = x - sup(j,:)';
        
        % RLS update
        Pr = P * r;
        k  = Pr / (1 + r' * Pr);
        P  = P - k * (r' * P);
        % Update readout weights
        phi = phi - k*e'; 
      end
    end
end
end


end


end