Solve system of equations (nleqslv) only returns origin

I'm solving 4 complex equations simultaneously.  Code is below.  The code
returns only zero's for the solution though there should also be a non-zero
result.  I'm pretty confident that the equations are correct because they
are straight from a published paper and I checked them pretty thoroughly.
The parameter values I used are from the published paper as well.  Any
suggestions for how I can find the maximum non-zero solution instead of
going straight to the minimum?

Thanks!
Alicia

install.packages("nleqslv",
lib="Users/Alicia/Documents/R.Codes/R.Packages/")

require(nleqslv)

###### Global Parameters ############
beeta=0.8
pq=10000
L=12600
theta=0.6
psale=0.6
mu=psale*(1-theta)
alphah=0.15
Cg=6240
Cs=2820
A= 100
D=0.0001
greekp=0.43
K=100000

##### Species Parameters ##########

b1=0.38
p1=16654
v1 = 0.28
N1=6000
g1=1
delta1=1

b2=0.4
p2=2797
v2 = 0.31
N2=10000
g2=1
delta2=1

### Define functions with vector x = c(Lg, Ls, gamma1, gamma2, lamda)

firstordercond <- function (x) {

   y=numeric(4)


y[1]=(alphah/x[3])-(x[5]*((p1-(((theta+mu)*(((N1/A)*g1^greekp*x[1]^b1)+K))+((theta+mu)*(((1-exp(-2*D*v1*N1))*x[2])+K))))*(((N1/A)*g1^(greekp))*x[1]^b1+(2*v1*N1*D)*x[2])
+
     delta1*theta*(((N1/A)*g1^(greekp))*x[1]^b1+(2*v1*N1*D)*x[2])))


y[2]=(alphah/x[4])-(x[5]*((p2-(((theta+mu)*(((N2/A)*g2^greekp*x[1]^b2)+K))+((theta+mu)*(((1-exp(-2*D*v2*N2))*x[2])+K))))*(((N2/A)*g2^(greekp))*x[1]^b2+(2*v2*N2*D)*x[2])
+
     delta2*theta*(((N2/A)*g2^(greekp))*x[1]^b2+(2*v2*N2*D)*x[2])))


y[3]=((alphah*((N1/A)*g1^(greekp))*b1*x[1]^(b1-1))/(((N1/A)*g1^(greekp))*x[1]^b1+(2*v1*N1*D)*x[2]))
+
((alphah*((N2/A)*g2^(greekp))*b2*x[1]^(b2-1))/(((N2/A)*g2^(greekp))*x[1]^b2+(2*v2*N2*D)*x[2]))
+
     x[5]*(((-1*beeta*pq*(L-x[1]-x[2])^(beeta-1)) +
(b1*(1-x[3])*(p1-(((theta+mu)*(((N1/A)*g1^greekp*x[1]^b1)+K))+((theta+mu)*(((1-exp(-2*D*v1*N1))*x[2])+K))))*((N1/A)*g1^(greekp))*x[1]^(b1-1))
+

(b2*(1-x[4])*(p2-(((theta+mu)*(((N2/A)*g2^greekp*x[1]^b2)+K))+((theta+mu)*(((1-exp(-2*D*v2*N2))*x[2])+K))))*((N2/A)*g2^(greekp))*x[1]^(b2-1))
-
   (delta1*theta*x[3]*((N1/A)*g1^(greekp))*b1*x[1]^(b1-1)) -
(delta2*theta*x[4]*((N2/A)*g2^(greekp))*b2*x[1]^(b2-1)))-Cg*(((N1/A)*g1^(greekp))*b1*x[1]^(b1-1)+((N2/A)*g2^(greekp))*b2*x[1]^(b2-1)))


y[4]=((alphah*(2*v1*N1*D))/(((N1/A)*g1^(greekp))*x[1]^b1+(2*v1*N1*D)*x[2]))
+ ((alphah*(2*v2*N2*D))/(((N2/A)*g2^(greekp))*x[1]^b2+(2*v2*N2*D)*x[2])) +
     x[5]*(((-1*beeta*pq*(L-x[1]-x[2])^(beeta-1)) +
((1-x[3])*(p1-(((theta+mu)*(((N1/A)*g1^greekp*x[1]^b1)+K))+((theta+mu)*(((1-exp(-2*D*v1*N1))*x[2])+K))))*(2*v1*N1*D))
+

((1-x[4])*(p2-(((theta+mu)*(((N2/A)*g2^greekp*x[1]^b2)+K))+((theta+mu)*(((1-exp(-2*D*v2*N2))*x[2])+K))))*(2*v2*N2*D))
- (delta1*theta*x[3]*(2*v1*N1*D)) -
     (delta2*theta*x[4]*(2*v2*N2*D)))-Cs*((2*v1*N1*D)+(2*v2*N2*D)))

     result=(x)

   }

Xstart=c(10000, 200, 0.5, 0.5, 12)
fstart= firstordercond(Xstart)

nleqslv(Xstart, firstordercond)



-- 
Alicia Ellis
Postdoc
Gund Institute for Ecological Economics
University of Vermont
617 Main Street
Burlington, VT  05405
(802) 656-1046
http://www.uvm.edu/~aellis5/
<http://entomology.ucdavis.edu/faculty/scott/aellis/>

	[[alternative HTML version deleted]]

You should return the vector y i.e. function values.
But y has length 4 and x has length 4.

So where is the fifth value for y?

  }

Xstart=c(10000, 200, 0.5, 0.5, 12)
fstart= firstordercond(Xstart)

If you print fstart you will see that it is identical to Xstart.

You need to rethink you firstordercond() function.
It's not correct.

Berend

nleqslv(Xstart, firstordercond)



-- 
Alicia Ellis
Postdoc
Gund Institute for Ecological Economics
University of Vermont
617 Main Street
Burlington, VT  05405
(802) 656-1046
http://www.uvm.edu/~aellis5/
<http://entomology.ucdavis.edu/faculty/scott/aellis/>

	[[alternative HTML version deleted]]