Re: steady state finder! still works after 6 hours!

[Nir Krakauer]

Do you even need symbolic math if all you want is a numeric solution? The system can be solved numerically in Octave with no extra packages:

pbeta=.99;

epsilon=-1; 

sigma=.03; 

theta=.3;

k=13;

delta=.025;

rho=1.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]);

#v = 1.1556   146.1857