For single-valued complex functions of a single complex variable
Most math software (free to use or proprietary) with 3D capabilities will have some form of built-in complex number arithmetic. While commands like plot(z,f(z)) may not work for complex $z$, it is usually possible to plot the real or imaginary parts of a complex-valued function by splitting the input into its real and imaginary components, and inputting <plot / surface>(x,y,real(f(x+iy))) and <plot / surface>(x,y,imag(f(x+iy))), respectively, or equivalent in place of $f(z)$.
In GeoGebra, for example, you can plot the real part of complex function $f$ by using $g(x,y)=\text{real}(f(x+iy))$. The exact syntax for this sort of thing will vary depending on the program you are using, but most programs will have a command to extract the real or imaginary part from a complex number.
Programs for which this works include MatLab, SciLab, Maple, Mathematica, GeoGebra, and anything that is built off of Python.
For multi-valued complex functions of a single complex variable
I've found that multi-valued functions are best realized as parametric surfaces. For a multivalued function $f$, you can start to find a parameterization by taking the inverse of $f$ - in practice, this amounts to writing $f(f^{-1}(z))=z$ (substituting the formula for $f$) and solving for $f^{-1}(z)$. You can then write the parametric surface as $\langle\Re[f^{-1}(u+iv)],\Im[f^{-1}(u+iv)],u\rangle$ for the real part and $\langle\Re[f^{-1}(u+iv)],\Im[f^{-1}(u+iv)],v\rangle$ for th imaginary.
As an example, suppose you want the Riemann surface for $f(z)=\sqrt{z}$. Solving for $f^{-1}$, you get
$f(z)=\sqrt{z}\to z=\sqrt{f^{-1}(z)}\to z^2=f^{-1}(z)$
$f^{-1}(u+iv)=(u+iv)^2=u^2+2uv\cdot i-v^2$
Riemann surface for the real part: $\langle u^2-v^2,2uv,u\rangle$
Programs for which this works include GeoGebra, Mathematica, Maple, and Math3D. While it's still possible in more technically oriented programs (like MatLab and SciLab), it does require a bit of coding to get it to work.