Unable to Solve steady state for three agents

var Y K C r y1 k1 c1 n1 a1 pi1 y2 k2 c2 n2 a2 pi2 y3 k3 c3 n3 a3 pi3 ek1 ek2 ek3;                                             
varexo eta_k1 eta_k2 eta_k3  ;
parameters theta2  Phi1 Phi2 Phi3 rho2 a0 pi0   
wy1 wy2 wy3 wc1 wc2 wc3 wk1 wk2 wk3 yc kc beta delta Psi1 Psi2 alpha phi2 phi3 
lambda1 lambda2 lambda31 lambda32 lambda33      
phi1 rho1 rho_k1 rho_k2 rho_k3;   

theta2 = 0.4970;   
Phi1 = 0.4930;     
Phi2 = 0.5670;     
Phi3 = 0.4870;     
rho2 = 0.0023;    

wy1 = 0.3;                    
wy3 = 0.2;                      
wy2 = 1-wy1-wy3;               
wk1 = 0.5;                   
wk3 = 0.2;                     
wk2 = 1-wk1-wk3;               
wc1 = 0.4;                     
wc3 = 0.2;                   
wc2 = 1 - wc1-wc3;             

yc = 2.3;                    % steady state of Y/C
kc = 16;                    % steady state of K/C
yk = yc/kc;                    % Steady state of Y/K
a0 = 0.016;                     
pi0 = 0.369;                     
 
alpha  = 0.3;                  
beta = 0.97;                    
delta = 0.06;                  
rss = -log(beta);          
theta1 = 0.016;              
Psi1 = 1;                     
Psi2 = 1;                      

phi1 = 0.9;                               
phi2 = alpha * yk / (delta + rss);        
phi3 = 1/(delta + rss) - phi2/Psi1;       
lambda1 = 1/(1+Psi2);                       
lambda2 = Psi1/(1+Psi2);                    
lambda31 = 2*Psi2*Phi1/(theta2*(1+Psi2)); 
lambda32 = 2*Psi2*Phi2/(theta2*(1+Psi2));             
lambda33 = 2*Psi2*Phi3/(theta2*(1+Psi2));  


rho1 = 0.9;             
rho_k1=0.95;           
rho_k2=0.95;             
rho_k3=0.95;         
sdeta_k1=0.95;      
sdeta_k2=0.95;        
sdeta_k3=0.95;          

model; 
	Y = wy1 * y1 + wy2 * y2 + wy3 * y3 ;   
    K = wk1 * k1 + wk2 * k2 + wk3 * k3;      l
    C = yc * Y - kc * (K-(1-delta)*K(-1)) ;   
    r = Psi1*(c1(+1)-c1)-log(beta);

    y1 = alpha * k1 + (1 - alpha) * n1+ a1;                            
    k1 = phi1 * k1(-1)+ phi2 * y1 - phi3 * r +ek1;                               
    c1= 1/wc1*(C -wc2*c2-wc3*c3)  ;                                   
    n1 = lambda1 * y1 - lambda2 * c1 + lambda31 * pi1 ;                     
    a1 = a0 + a1(-1) - Phi1 * pi1(-1) ;                                       
    pi1 = pi0 + rho1 * pi1(-1) - rho2 * k1(-1);                             
    ek1 = rho_k1 * ek1(-1)+eta_k1;


    y2 = alpha * k2 + (1 - alpha) * n2 + a2 ;                          
    k2 = phi1 * k2(-1)+ phi2 * y2 - phi3 * r + ek2;                  
    c2=c2(+1)-Psi1*(r+log(beta))   ;                         
    n2 = lambda1* y2 - lambda2 * c2 + lambda32 * pi2;                  
    a2 = a0 + a2(-1) +Phi2 * pi2(-1) ;                                        
    pi2 = pi0 + rho1 * pi2(-1) - rho2 * k2(-1) ;                             
    ek2 = rho_k2 * ek2(-1)+eta_k2;

    y3 = alpha * k3 + (1 - alpha) * n3 + a3;                            
    k3 = phi1 * (y3 + 1/Psi1 * log(beta)) + phi3 * r+ek3;                    
    c3=c3(+1)-Psi1*(r+log(beta))    ;                                        
    n3 = lambda1* y3 - lambda2 * c3 + lambda33 * pi3;                     
    a3 = a0 + a3(-1) + Phi3 * pi3(-1) ;                                       
    pi3 = pi0 + rho1 * pi3(-1) - rho2 * k3(-1) ;                              
    ek3 = rho_k3 * ek3(-1)+eta_k3;
end; 

initval;
 Y=0;
 K=0 ;
 C=0 ;
 r=rss;
 y1=0; 
 k1=0; 
 c1=0;
 n1=0; 
 a1=0; 
 pi1=0; 
 y2=0; 
 k2=0; 
 c2=0; 
 n2=0;
 a2=0; 
 pi2=0; 
 y3=0; 
 k3=0; 
 c3=0; 
 n3=0; 
 a3=0; 
 pi3=0; 
 ek1=0 ;
 eta_k1=0.25 ;
 ek2=0 ;
 eta_k2=0.25 ; 
 ek3=0 ;
 eta_k3=0.25 ;
end; 

shocks;
var eta_k1;
stderr sdeta_k1;
var eta_k2;
stderr sdeta_k2;
var eta_k3;
stderr sdeta_k3;
end;

steady;
resid;
%stoch_simul(periods = 200,irf = 40,order=1);
%dynasave('simudata.mat');
%figure(1)
%subplot(3,3,1)
% plot(oo_.endo_simul(1,:))

Your model is linearized, but some equations contain constant terms (like the one for r). You either need to leave out all constants or consistently include them.

Thanks a lot. I removed the constants as the following model block, it gives the steady states, but when I do the simulation, it gives the error message.

Blanchard Kahn conditions are not satisfied: no stable equilibrium

I use the check command, it gives

*There are 3 eigenvalue(s) larger than 1 in modulus *
for 1 forward-looking variable(s)
The rank condition ISN’T verified!

Could you please give me some suggestions?

model(linear); 
	Y = wy1 * y1 + wy2 * y2 + wy3 * y3 ;      % eq.1 Aggregate output
    K = wk1 * k1 + wk2 * k2 + wk3 * k3;      % eq.2 Aggregate capital
    C = yc * Y - kc * (k1-(1-delta)*k1(-1)) ;   % eq.3 Aggregate consuption = Output - Investment
    r = c2 (+1)- c2;       

    y1 = alpha * k1 + (1 - alpha) * n1+ a1;                            % eq.5 Rich's output
    k1= phi1 * k1(-1)+ phi2 * y1 - phi3 * r + ek1;                                % eq.6 Rich's capital
    c1= 1/wc1*(C -wc2*c2-wc3*c3)  ;                                     % eq.7 Consumption Euler equation of Rich
    n1 = lambda1 * y1 - lambda2 * c1 + lambda31 * pi1 ;                      % eq.8 Rich's labour supply
    a1 = a1(-1) - Phi1 * pi1(-1);                                       % eq.9 Rich's productivity
    pi1 =  rho1 * pi1(-1) - rho2 * k1(-1);                               % eq.10 Rich's penalty rate
    ek1 = rho_k1 * ek1(-1)+eta_k1;

    y2 = alpha * k2 + (1 - alpha) * n2 + a2;                            
    k2 = phi1 * k2(-1)+ phi2 * y2 - phi3 * r + ek2;                     
    c2 = yc * y2 - kc * (k2-(1-delta)*k2(-1));                          
    n2 = lambda1* y2 - lambda2 * c2 + lambda32 * pi2;                  
    a2 = a2(-1) - Phi1 * pi2(-1);                                      
    pi2 = rho1 * pi2(-1) - rho2 * k2(-1) ;                            
    ek2 = rho_k2 * ek2(-1)+eta_k2;

    y3 = alpha * k3 + (1 - alpha) * n3 + a3;                             
    k3 = phi1 * k3(-1)+ phi2 * y3 - phi3 * r +ek3;                      
    c3 = yc * y3 - kc * (k3-(1-delta)*k3(-1))   ;                          
    n3 = lambda1* y3 - lambda2 * c3 + lambda33 * pi3;                  
    a3 = a3(-1) + Phi3 * pi3(-1) ;                                       
    pi3 = rho1 * pi3(-1) - rho2 * k3(-1) ;                         
    ek3 = rho_k3 * ek3(-1)+eta_k3;
end;

Are you sure the timing of capital is correct, e.g. in the production function?