function [ F ]  = FindSS(x,alfa,omega,delta,csi,rkB,rkG,pI,G,KG,pl,taubar,ps,zeta,deltax,erow,D0,D1,D2)
y=x(1);
pB=x(2);
k=x(3);
e=x(4);

kG=KG*k;
kB=k-kG;
    yB = (rkB * kB) / (alfa * pB);
    for i = 1:1000  % Iterative approximation (Newton’s method could be used here)
        yB_new = (rkB * kB) / (alfa * (pB - taubar *(1-ps)* (yB^(1-ps))));
        if abs(yB_new - yB) < 1e-20  % Convergence condition
            break;
        end
        yB = yB_new;
    end

X=(e+erow)/deltax;
A=1-D0-D1*X-D2*X^2;
pG=(1/(1-zeta)*(pI^(1-csi)-zeta*((pB)^(1-csi))))^(1/(1-csi));
yG=(1-zeta)*(pG/pI)^(-csi)*y;   
Lb=(yB/(A*kB^(alfa)))^(1/(1-alfa));
Lg=(yG/(A*kG^(alfa)))^(1/(1-alfa));
iB=delta*kG;
iG=delta*kB;
i=iG+iB;                           
wG=(1-alfa)*pG*yG/Lg;
L=(Lb^(1+pl)+Lg^(1+pl))^(1/(1+pl));
g=G*y;

c=y-i-g; 

wB=(1-omega)/omega*c/((1-L)*(Lb^(1+pl)+Lg^(1+pl))^(pl/(1+pl))*Lb^(-pl));

F=[wB - (1 - alfa) * yB * (pB - taubar *(1-ps)* yB^(1 - ps)) / Lb 
    (1-L)*(Lb^(1+pl)+Lg^(1+pl))^(pl/(1+pl))*Lg^(-pl)- (1-omega)/omega*c/wG
    rkG*kG-alfa*pG*yG
    e-yB^(1-ps)];

end


