title: Sci Py Finding a set of nonnegative roots to system of equations
body: <p>I have a system of non-linear equations. I do not have a good guess about initial values. And I want at least one set of all positive roots cause in economics, negative values to these variables do not make much sense.   </p>
  
  <pre><code># -*- coding: utf-8 -*-
  """
  Created on Sat Oct 15 21:48:56 2016
  
  @author: Nick
  """
  
  import scipy as sp
  from scipy.optimize import root, fsolve
  import numpy as np
  
  #from scipy.optimize import *
  
  
  
  el          = 1.1  
  eg          = el   
  ej          = 10  
  om          = 0.3  
  omg         = 0.3  
  rhog        = 0.8  
  xi          = 0.9  
  mun         = 2
  pidss       = 0.02 
  muc         = 0.001
  ec          = 2.00 # sims obtains 2.47
  beta        = 0.998
  h           = 0.8   
  kappa       = 4.00 
  n           = 1/3.0  
  alpha       = 1/3.0  
  delta       = 0.025
  egs         = eg   
  oms         = 0.2  
  omgs        = oms  
  rhom        = 0.7  
  psiygap     = 1.000
  psipi       = 2.500
  rhoicu      = 0.800
  taudss      = 0.01    # steady state tax on domestic consumption (setting it as 0 would create algebraic difficulties)
  taumss      = 0.01    # steady state tax on imported consumption for domestic country
  taukss      = 0.01    # steady state tax on rental income from capital for domestic country block
  taunss      = 0.01    # steady state tax on labor for domestic country
  tauydss     = 0.05    
  gss         = 0.23    # steady state government spending as a propostion of gdp for domestic country block    
  gsss        = 0.23    # steady state government spending as a propostion of gdp for foreign country block    
  
  taudsss     = 0.01    
  
  taumsss     = 0.01    
  
  tauksss     = 0.01    
  
  taunsss     = 0.01    
  
  tauydsss    = 0.01          # steady state tax rate on output for foreign country block
      tauss       = 1.0               # Steady state terms of trade
  
  icu = ((1+pidss)/beta) - 1 
  mc = ((ej - 1)/ej)
  r = (1/taukss) * ((1/beta)  - (1-delta))
  rs = (1-tauksss) * ((1/beta)  - (1-delta))
  KN = (mc*alpha/r)**(1/(1-alpha))
  KNs = (mc*alpha/rs)**(1/(1-alpha))
  psigma = (1-xi) * (1/(1-tauydss) - xi)**(-1)
  psigmas = (1-xi) * (1/(1-tauydsss) - xi)**(-1)
  w =     (1-alpha) * mc * (KN)**(alpha)
  
  z = np.zeros(16)
  
  def fun(z):
      Yd = z[0]
      N = z[1]
      X = z[2]
      I = z[3]
      Cd = z[4]
      Cm = z[5]
      Gd = z[6]
      Gm = z[7]
      Yds = z[8]
      Ns = z[9]
      Xs = z[10]
      Is = z[11]
      Cds = z[12]
      Cms = z[13]
      Gds = z[14]
      Gms = z[15]
      print (z)
      f = np.zeros(16)
      f[0]  = N - ( (X - muc)**(-ec) * ((1-alpha)/(mun)) * (mc)**(1/(1-alpha)) * (alpha/r)** (1-taunss) )
      f[1]  = Yd - ( Cd + Gd + I + ((1-n)/n) *(Cms + Gms)       )
      f[2]  = Yd - ( (KN)**(alpha)  * (psigma/(1-tauydss)**(ej)) )
      f[3]  = Cd - ( X * ((1-om)/(1+taudss)**(el)) *((1-om)*(1+taudss)**(1-el) + om * (1+taumss)**(1-el) * tauss**(1-el))**(el/(1-el))  )
      f[4]  = Gd - ( ((gss*Yd * (1-omg))/(1+taudss)**(eg) ) *((1-omg)*(1+taudss)**(1-eg) + omg* (1+taumss)**(1-eg) * tauss**(1-eg))**(eg/(1-eg) ) )
      f[5]  = I - ( delta* KN * N )
      f[6]  = Cm -( (X * (1-om)/(1+tauydss)**(el) ) *((1-om)*(1+taudss)**(1-el) + om* (1+taumss)**(1-el) * tauss**(1-el))**(el/(1-el))   )
      f[7]  = Gm - ( ((gss*Yd * (omg))/(1+taumss)**(eg) ) *((1-omg)*(1+taudss)**(1-eg) + omg* (1+taumss)**(1-eg) * tauss**(1-eg))**(eg/(1-eg) ) )
      f[8]  = Ns - ( (Xs - muc)**(-ec) * ((1-alpha)/(mun)) * (mc)**(1/(1-alpha)) * (alpha/rs)** (1-taunsss) )
      f[9]  = Yds - ( Cds + Gds + Is + (n/(1-n)) *(Cm + Gm)       )
      f[10] = Yds - ( (KNs)**(alpha)  * (psigmas/(1-tauydsss)**(ej) ) )
      f[11] = Cds - ( Xs * ((1-oms)/(1+taudsss)**(el))* ((1-oms)*(1+taudsss)**(1-el) + oms* (1+taumsss)**(1-el) * tauss**(1-el))**(el/(1-el))  )
      f[12] = Gds - ( ((gsss*Yds * (1-omgs))/(1+taudsss)**(eg) ) *((1-omgs)*(1+taudsss)**(1-eg) + omgs* (1+taumsss)**(1-eg) * tauss**(1-eg))**(eg/(1-eg) ) )
      f[13] = Is - ( delta* KNs * Ns )
      f[14] = Cms -( (Xs * (1-oms)/(1+tauydsss)**(el) ) *((1-oms)*(1+taudsss)**(1-el) + oms* (1+taumsss)**(1-el) * tauss**(1-el))**(el/(1-el))   )
      f[15] = Gms - ( ((gsss*Yds * (omgs))/(1+taumsss)**(eg) ) *((1-omgs)*(1+taudsss)**(1-eg) + omgs* (1+taumsss)**(1-eg) * tauss**(1-eg))**(eg/(1-eg) ) )
      return f
  
  
  
  z = sp.optimize.root(fun, [100,100,70,30,50,20,50,20,100,100,100,100,100,100,100,100],  method='lm')
  #z = fsolve(fun, [0,0,0.0,0,1,1,1,1,1,1,1,1,1,1,1,1])    
  print(z)
  </code></pre>
  
  <p>The solution are the following roots</p>
  
  <pre><code>success: True
         x: array([  3.64725445e-01,   1.02848541e-06,  -1.86761721e+02,
           9.52089296e-10,  -1.30733205e+02,  -1.25265418e+02,
           5.87207967e-02,   2.51660557e-02,   3.36422990e+00,
           5.18324506e-04,   8.17060628e+01,   4.87111630e-04,
           6.53648502e+01,   6.53648502e+01,   6.19018302e-01,
           1.54754576e-01])
  </code></pre>