title: Python: multivariate non-linear solver with constraints
body: <p>Given a function <code>f(x)</code> that takes an input vector <code>x</code> and returns a vector of the same length, how can you find the roots of the function setting constraints on <code>x</code>? (E.g. a range for each component of <code>x</code>.)</p>
  
  <p>To my surprise I could not find a lot of useful information about this. In the scipy list for <a href="https://docs.scipy.org/doc/scipy/reference/optimize.html#module-scipy.optimize" rel="noreferrer">Optimization and Root finding algorithms</a> there seem to be some options for scalar functions such as <a href="https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.brentq.html#scipy.optimize.brentq" rel="noreferrer">brentq</a>. I can not find any algorithms that supports such an option for the multivariate case though.</p>
  
  <p>Of course one could do a work-around like squaring each component of the returned vector and then use one of the minimizers such as <a href="https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.differential_evolution.html#scipy.optimize.differential_evolution" rel="noreferrer">differential_evolution</a> (this is the only one I think actually). I can not imagine that this is a good strategy though, since it kills the quadratic convergence of Newton's algorithm. Also I find it really surprising that there does not seem to be an option for this, since it must be a really common problem. Have I missed something?</p>