1 | initial version |

You could use mpmath's multidimensional solver, but unfortunately it's not well-wrapped either. On the bright side you can do it without any Cython:

```
import mpmath
var("a b")
eq0 = -2*(a*b/sqrt(-1/2*a^2 + b^2) - 2/(a^2*b^2))*(sqrt(-1/2*a^2 + b^2)*b - 1/(a*b^2)) - 16/(a^5*b^2)==0
eq1 = 4*b^3 + 4*(sqrt(-1/2*a^2 + b^2)*b - 1/(a*b^2))*(b^2/sqrt(-1/2*a^2 + b^2) + sqrt(-1/2*a^2 + b^2) + 2/(a*b^3)) - 8/(a^4*b^3)==0
f = [lambda a,b: eq0.lhs().subs(a=RR(a), b=RR(b)),
lambda a,b: eq1.lhs().subs(a=RR(a), b=RR(b))]
found_root = mpmath.findroot(f, (2, 2))
found_root = Matrix(RR, found_root.tolist())
print found_root
fa,fb = found_root.list()
print eq0.subs(a=fa,b=fb)
print eq1.subs(a=fa,b=fb)
```

which gives:

```
sage: load "rfind.sage"
[1.49174202118713]
[1.10177415033843]
(-2.67785793539588e-13) == 0
(4.93827201353270e-13) == 0
```

and so we see it's done a decent job. The above was a very manual wrapping, it wouldn't be that much harder to treat the general case.

Copyright Sage, 2010. Some rights reserved under creative commons license. Content on this site is licensed under a Creative Commons Attribution Share Alike 3.0 license.