SymPy Problem Instance

import sympy
from sympy.solvers.solveset import nonlinsolve
dz = 4.132000000e+16
x0 = [
    0.0          , #theta
    1064 * 1.0e-9, #lambda
    5.0          , #r_w
    1.05         , #M2
    0.0          , #f_a
    0.0           #r_s
]
theta, lamb, r_w, M2, f_a, r_s = sympy.symbols("a, b, c, d, x, y", real=True)
sol = nonlinsolve([
    theta - M2*lamb/(sympy.pi*r_w),
    r_w*r_w/((r_w + dz*sympy.tan(theta))**2) - f_a,
    (r_w + dz*sympy.tan(theta)) - r_s
], [theta, lamb, r_w, M2, f_a, r_s])