For a gaussian beam, how can I calculate $w_0$, $z$ and $z_R$ from knowledge of $w$, $R$ and $\lambda$? I have a Gaussian beam somewhere away from the beam waist (also referred to as focus point). Then, according to text books (or wikipedia), the beam size $w$ is given by 
$$
w(z) = w_0 \sqrt{ 1+\left(\frac{z}{z_R}\right)^2 },
$$
with $z$ the axial distance from the beam waist, $w_0$ the waist size and the Rayleigh length $z_R$ which is defined as 
$$
z_R=\frac{\pi w_0^2}{\lambda_0}.
$$
The radius of curvature of the phase fronts are given by
$$
R(z) = z \left( 1 + \left(\frac{z_R}{z}\right)^2 \right).
$$
I have: $R$, $w$, and $\lambda_0$.
I am looking for: $w_0$, $z$, and $z_R$. 
My idea was to solve the third equation for $z$,
$$
z_{1,2} = \frac{R}{2} \pm \sqrt{\left( \frac{R}{2} \right)^2 - z_R^2},
$$
and then inserting into the first equation 
$$
w = w_0 \sqrt{ 1 + \left( \frac{ \frac{R}{2} \pm \sqrt{\left( \frac{R}{2} \right)^2 - z_R^2} }{z_R} \right)^2 }.
$$
Inserting $z_R$ yields
$$
w = w_0 \sqrt{ 1 + \left( \frac{ \frac{R}{2} \pm \sqrt{\left( \frac{R}{2} \right)^2 - \left(\frac{\pi w_0^2}{\lambda_0}\right)^2} }{\frac{\pi w_0^2}{\lambda_0}} \right)^2 }
$$
which is an equation for $w_0$. Since I am currently stuck solving it, my question is in principle twofold: what software/tool do you guys recommend me for solving such an equation for $w_0$ and if there is a different way determining the beam parameters (maybe I missed something) ?
 A: When attacking a system of equations with square roots it is often helpful to get rid of the square roots until the very end (if possible). Good idea garyp. 
Let's start with...
$$w^2 = w_0^2 \left( 1 + \frac{z^2}{z_R^2} \right)$$
$$z_R = \frac{\pi w_0^2}{\lambda_0}$$
$$R = z \left( 1 + \frac{z_R^2}{z^2}\right)$$
When I look at this system of equations I notice that $z_R = \pi w_0^2 / \lambda_0$ is only relating $z_R$ and $w_0$, two unknowns, with $\lambda_0$, a known.  This means I can easily use $z_R = \pi w_0^2 / \lambda_0$ to replace either the $w_0$ or the $z_R$ in all of the equations and I would be down to two equations with two unknowns. I recommend getting rid of the $w_0^2$ first as it only shows up in one place.* 
$$w^2 = \frac{\lambda_0 z_R}{\pi}\left( 1 + \frac{z^2}{z_R^2} \right)$$
$$R = z \left( 1 + \frac{z_R^2}{z^2}\right)$$
Now move the $\lambda_0/\pi$ to the left hand side and expand the right hand sides to see that there is symmetry in the two remaining equations. 
$$ \frac{\pi w^2}{\lambda_0 }  =z_R  + \frac{z^2}{z_R}$$
$$R = z  + \frac{z_R^2}{z}$$
Notice that if I make both right hand sides the same we can set them equal to each other. 
$$ \frac{\pi w^2}{\lambda_0 }z_R  =z_R^2  + z^2 $$
$$R z = z^2  + z_R^2$$
Now it is easy to use $ z =  (\pi w^2/(\lambda_0 R))z_R$ to plug into the  earlier forms . 
$$ \frac{\pi w^2}{\lambda_0 }  =z_R  + \frac{z^2}{z_R} = z_R \left( 1  + \frac{\pi^2 w^4}{\lambda_0^2 R^2} \right)  $$
$$R = z  + \frac{z_R^2}{z} = z \left( 1 + \frac{ \lambda_0^2 R^2}{\pi^2w^4} \right) $$
The final results are 
$$ z_R  = \frac{ \pi \lambda_0 R^2 w^2 } {\lambda_0^2 R^2  + \pi^2 w^4}   $$
$$ z =\frac{\pi^2 R  w^4}{ \lambda_0^2 R^2 + \pi^2 w^4} $$
and we can use the relation for $z_R$ and $w_0$ to get 
$$  w_0 = \frac{  \lambda_0 \left|R \right| w } {\sqrt{\lambda_0^2 R^2  + \pi^2 w^4}} $$
Now $w_0$ must be greater than $0$ and so must $w$ but $R$ could be positive or negative. Therefore we need an absolute value on the $R$ in the above equation.  
heather was right in that Mathematica would be very helpful. 
Solve[{w^2 == w0^2 (1 + z^2/zr^2), zr == \[Pi] w0^2/lambda0, R == z (1 + zr^2/z^2)}, {w0, z, zr}] // FullSimplify

*I would like to mention that this step wasn't obvious to me from the start. I figured it out by realizing that $\pi/\lambda_0$ was not that important so I set $\lambda=\pi$. This way  I could see the key features of the equations without the mess of coefficients.
A: The easiest way to do this is via the complex beam parameter $q$, which is given by
$$
\frac1q = \frac{1}{z+iz_R} = \frac{1}{R(z)}-i\frac{\lambda_0}{\pi w(z)^2}.
$$
In your case you know all three parameters on the right hand side, so you know $q$, and therefore you can determine $z$ and $z_R$ as its real and imaginary parts. Doing this explicitly,
\begin{align}
q
& =
\frac{1}{\frac{1}{R(z)}-i\frac{\lambda_0}{\pi w(z)^2}}
\\ & =
\frac{ R(z)\pi w(z)^2}{\pi w(z)^2-i\lambda_0R(z)}
\\ & =
\frac{ R(z)\pi w(z)^2}{\pi w(z)^2-i\lambda_0R(z)} \frac{\pi w(z)^2+i\lambda_0R(z)}{\pi w(z)^2+i\lambda_0R(z)}
\\ & =
\frac{ R(z)\pi w(z)^2}{\pi^2 w(z)^4+\lambda_0^2R(z)^2}\left(\pi w(z)^2+i\lambda_0R(z)\right),
\end{align}
and from there you can read off
$$
z=
\frac{R(z)\pi^2 w(z)^4}{\pi^2 w(z)^4+\lambda_0^2R(z)^2}
$$
and
$$
z_R=
\frac{\lambda_0R(z)^2\pi w(z)^2}{\pi^2 w(z)^4+\lambda_0^2R(z)^2}.
$$
Finally, the beam waist is determined by the Rayleigh range:
$$
w_0^2
=\frac{\lambda_0}{\pi}z_R 
= \frac{\lambda_0^2 R(z)^2 w(z)^2}{\pi^2 w(z)^4+\lambda_0^2R(z)^2},
$$
so
$$
w_0
= \frac{\lambda_0 R(z) w(z)}{\sqrt{\pi^2 w(z)^4+\lambda_0^2R(z)^2}}.
$$
These agree with LasersMatter's results, but the $q$ route is a lot easier.
