I have the following system of non-linear equations and want to find their solutions experimenting with my parameters. I also want to restrict the solutions to be non-negative. I have done the following, but i am sure it exist a more efficient way. Can somone help on this?

eqns := [A = (gr_c+delta)*kh^(1-alpha)/sav_rate, theta = Rk*(Rh-Rk)/(gamma*((Rh-Rk)^2+sigma^2)), theta = 1*1+kh, Rk = 1+rk-delta, Rh = 1+rh-delta, rk = A*alpha*kh^(alpha-1), rh = A*(1-alpha)*kh^alpha, sigma = sigmay/theta, varrho = Rap^((eps-1)*eps/(1-gamma)), Rap = Rk^(1-gamma)+(1-gamma)*Rk^(-gamma)*theta*(Rh-Rk)-.5*Rk^(-gamma-1)*gamma*(1-gamma)*theta^2*((Rh-Rk)^2+sigma^2), R = Rk+theta*(Rh-Rk), beta = ((1+gr_c)/R)^(1/eps)/varrho];

print(`output redirected...`); # input placeholder

[

[

[

[

[ (1 - alpha)

[ (gr_c + delta) kh

[A = ----------------------------,

[ sav_rate

[

Rk (Rh - Rk)

theta = ---------------------------, theta = 1 + kh,

/ 2 2\

gamma \(Rh - Rk) + sigma /

Rk = 1 + rk - delta, Rh = 1 + rh - delta,

(alpha - 1) alpha

rk = A alpha kh , rh = A (1 - alpha) kh ,

/(eps - 1) eps\

|-------------|

sigmay \ 1 - gamma /

sigma = ------, varrho = Rap , Rap =

theta

(1 - gamma) (-gamma)

Rk + (1 - gamma) Rk theta (Rh - Rk) - 0.5

(-gamma - 1) 2 / 2 2\

Rk gamma (1 - gamma) theta \(Rh - Rk) + sigma /,

/ 1 \]

|---|]

\eps/]

/1 + gr_c\ ]

|--------| ]

\ R / ]

R = Rk + theta (Rh - Rk), beta = ---------------]

varrho ]

]

vals := [alpha = .36, delta = 0.6e-1, sigmay = sqrt(0.313e-1), gamma = 3, eps = .5, gr_c = 0.2e-1, sav_rate = .23];

eval(eqns, vals);

print(`output redirected...`); # input placeholder

[

[

[

[ 0.64 Rk (Rh - Rk)

[A = 0.3478260870 kh , theta = -----------------------,

[ / 2 2\

[ 3 \(Rh - Rk) + sigma /

0.36 A

theta = 1 + kh, Rk = 0.94 + rk, Rh = 0.94 + rh, rk = ------,

0.64

kh

0.36 0.1769180601

rh = 0.64 A kh , sigma = ------------,

theta

0.1250000000

varrho = Rap ,

1 2 theta (Rh - Rk)

Rap = --- - -----------------

2 3

Rk Rk

2 / 2 2\

3.0 theta \(Rh - Rk) + sigma /

+ --------------------------------, R = Rk + theta (Rh - Rk),

4

Rk

2.000000000]

/1\ ]

1.0404 |-| ]

\R/ ]

beta = ---------------------]

varrho ]

]

eqns := [A = .3478260870*kh^.64, theta = (1/3)*Rk*(Rh-Rk)/((Rh-Rk)^2+sigma^2), theta = 1+kh, Rk = .94+rk, Rh = .94+rh, rk = .36*A/kh^.64, rh = .64*A*kh^.36, sigma = .1769180601/theta, varrho = Rap^.1250000000, Rap = 1/Rk^2-2*theta*(Rh-Rk)/Rk^3+3.0*theta^2*((Rh-Rk)^2+sigma^2)/Rk^4, R = Rk+theta*(Rh-Rk), beta = 1.0404*(1/R)^2.000000000/varrho];

print(`output redirected...`); # input placeholder

[

[

[

[ 0.64 Rk (Rh - Rk)

[A = 0.3478260870 kh , theta = -----------------------,

[ / 2 2\

[ 3 \(Rh - Rk) + sigma /

0.36 A

theta = 1 + kh, Rk = 0.94 + rk, Rh = 0.94 + rh, rk = ------,

0.64

kh

0.36 0.1769180601

rh = 0.64 A kh , sigma = ------------,

theta

0.1250000000

varrho = Rap ,

1 2 theta (Rh - Rk)

Rap = --- - -----------------

2 3

Rk Rk

2 / 2 2\

3.0 theta \(Rh - Rk) + sigma /

+ --------------------------------, R = Rk + theta (Rh - Rk),

4

Rk

2.000000000]

/1\ ]

1.0404 |-| ]

\R/ ]

beta = ---------------------]

varrho ]

]

solve(eqns, [Rk, Rh, varrho, Rap, beta, R, A, sigma, theta, rk, rh, kh]);