function [ys,params,check] = example_PissaridesIV_steadystate(ys,exo,M_,options_)

% read out parameters to access them with their name
NumberOfParameters = M_.param_nbr;
for ii = 1:NumberOfParameters
  paramname = M_.param_names{ii};
  eval([ paramname ' = M_.params(' int2str(ii) ');']);
end
% initialize indicator
check = 0;

%% Enter model equations here

if length(exo) >= 1
    flow = exo(1);
end

F = F0 + flow;

r=(1 / bbeta) - 1;
rr =(1 - bbeta);

load('stst_guess.mat')
run solve_SS.m

f = mu * theta^(1 - eta);

U = (F * d)/(d + f);
L =  F - U;
v = U * theta;
q = f / theta;
Aj  = rr + bbeta * d;
Y   = (K^aalpha) * (L^(1 - aalpha)); % F(L,K)

J_b = (pMg - w) / Aj;
i=r;
I = ddd * K;
C = Y - (I + v * c);

syms W_b U_b

eqW_SS = W_b == (w * (1 - t) / C) + bbeta * ((1 - d) * W_b + d * U_b);
eqU_SS = U_b == (z / C) + bbeta * (f * W_b + (1 - f) * U_b);

eqsUW_SS=[eqW_SS, eqU_SS];

varsUW_SS=[W_b, U_b];
          [W_b, U_b]=vpasolve(eqsUW_SS, varsUW_SS);

params=NaN(NumberOfParameters,1);
for iter = 1:length(M_.params) %update parameters set in the file
  eval([ 'params(' num2str(iter) ') = ' M_.param_names{iter} ';' ])
end

NumberOfEndogenousVariables = M_.orig_endo_nbr; %auxiliary variables are set automatically
for ii = 1:NumberOfEndogenousVariables
  varname = M_.endo_names{ii};
  eval(['ys(' int2str(ii) ') = ' varname ';']);
end
 