function res = res_lss(l_grid, modelParams)

%parameters = [beta psi delta alpha rho sigma epsilon];
beta    = modelParams(1);
psi     = modelParams(2);
delta   = modelParams(3);
alpha   = modelParams(4);
rho     = modelParams(5);
sigma   = modelParams(6);
epsilon = modelParams(7);

r_ss = 1/beta-1+delta;
kl   = ( r_ss/((epsilon-1)/epsilon*alpha) )^(1/(alpha-1));
w_ss = (epsilon-1)/epsilon*(1-alpha)*kl^alpha;

%{
%plotting residuals of l for all the range (0,1):
l_grid = linspace(0,1,100);
res_l = (1-l_grid)*w_ss/psi + (delta*kl - kl^alpha)*l_grid;
plot(l_grid',[zeros(length(res_l),1) res_l'])
%}

res = (1-l_grid)*w_ss/psi + (delta*kl - kl^alpha)*l_grid;