function v = lucas4_steady_state_helper(pbeta)
options=optimset('Display','Final','TolX',1e-10,'TolFun',1e-10);
%addpath C:\dynare\4.4.3\matlab
gamma=2;
k=13;
epsilon=-1;
theta=.3;
sigma=.04;
delta=.01;
S=1;
rho=1;
psivar=.1;
Ybar=1;
Lambdabar=pbeta;
Rbar=1/pbeta;
Rebar=1/pbeta;
Rlbar= @(Qbar) 1 + 1 ./ Qbar;
munbar= @(Omegabar) Lambdabar*Omegabar*Rbar;
musbar= @(Omegabar, Qbar) Lambdabar*Omegabar*(Rlbar(Qbar)-Rbar);
muebar = @(Omegabar) Lambdabar*Omegabar*(1-pbeta);
epersbar = @(Omegabar, Qbar) musbar(Omegabar, Qbar) ./ muebar(Omegabar);
xbar = @(Omegabar, Qbar) -epersbar(Omegabar, Qbar)+(epersbar(Omegabar, Qbar)^(2)+(2/k)*(1-epsilon*epersbar(Omegabar, Qbar))).^(1/2);
Thetabar = @(Omegabar, Qbar) theta*(1+epsilon*xbar(Omegabar, Qbar)+(k/2)*xbar(Omegabar, Qbar).^(2));
phibar = @(Omegabar, Qbar) Omegabar ./ (Thetabar(Omegabar, Qbar)-musbar(Omegabar, Qbar)-xbar(Omegabar, Qbar).*muebar(Omegabar));
Ebar = @(Omegabar, Qbar) xbar(Omegabar, Qbar) .* Qbar;
Nbar = @(Omegabar,Qbar) Qbar ./ phibar(Omegabar, Qbar);
Dbar = @(Omegabar,Qbar) Qbar-Nbar(Omegabar,Qbar)-Ebar(Omegabar,Qbar);
Nnbar = @(Qbar) delta*(Qbar+Ybar);
eqs = @(Omegabar, Qbar) [Omegabar - (sigma+(1-sigma)*(munbar(Omegabar)+phibar(Omegabar, Qbar).*musbar(Omegabar, Qbar)+phibar(Omegabar, Qbar).*xbar(Omegabar, Qbar).*muebar(Omegabar))), Nbar(Omegabar,Qbar) - ((1-sigma).*(Qbar+Ybar-Rebar.*Ebar(Omegabar,Qbar)-Rbar*Dbar(Omegabar,Qbar))+Nnbar(Qbar))];
v = fsolve (@(x) eqs(x(1), x(2)), [1, 100]);