Orbitals of multi-electron atoms are qualitatively similar to those of hydrogen
Now that we have the closed-form analytical expressions for the Hydrogen wavefunctions, namely:
$${\displaystyle \psi _{n\ell m}(r,\theta ,\varphi )={\sqrt {{\left({\frac {2}{na_{0}^{*}}}\right)}^{3}{\frac {(n-\ell -1)!}{2n(n+\ell )!}}}}e^{-\rho /2}\rho ^{\ell }L_{n-\ell -1}^{2\ell +1}(\rho )Y_{\ell }^{m}(\theta ,\varphi )}$$
(where $\rho=2r/(na_\mathrm{B})$, $L$ are Laguerre polynomials $Y$ spherical harmonics) we can turn to the problem of their visualisation.
The orbitals represent the probability density function to find the electron according to the postulate of quantum mechanics. At such, they can be understood as a description of a diffuse electron cloud in a standing wave pinned to its nucleus.
They extend to much extent to all atoms, as we shall see in the coming lectures, where they are called "atomic orbitals". For hydrogen, the stationary orbitals are specified by three quantum numbers:
$$n,\quad l\quad\text{and}\quad m$$
The $s$ subshells are shaped like spheres, with the ground state being simply a fading ball. Higher $n$ with $l=0$ have radial (spherical) symmetry but with oscillations due to destructive interferences. The $s$ orbitals are the only ones with an antinode at the centre (where the nucleus is) while other orbitals present a node there. This does not invalidate the model, though, as the probability for the electron to be actually found there remain negligible.
There are three $p$-orbitals (starting at $n=2$) corresponding to $m=\pm1$ and 0. They each have the form of two ellipsoids pointing in opposite directions and pinned by a node at the nucleus, and are otherwise perpendicularly oriented the ones from the others, which one can align with the Cartesian axis vectors $\hat\imath$, $\hat\jmath$ and $\hat k$.
There are five $d$-orbitals (starting at $n=3$) corresponding to $m=\pm2, \pm1$ and 0.
(linear combination)
(see for instance this Wikipedia page). A great piece of code is provided by Jacopo Bertolotti:
\[Alpha]0 = 1; \[Psi][n_, l_, m_, r_, \[Theta]_, \[Phi]_] := Sqrt[(2/(n \[Alpha]0))^3 (n - l - 1)!/(2 n ((n + l)!))] E^(-r/(n \[Alpha]0)) ((2 r)/(n \[Alpha]0))^l LaguerreL[n - l - 1, 2 l + 1, (2 r)/(n \[Alpha]0)] SphericalHarmonicY[l, m, \[Theta], \[Phi]]; p1 = Flatten@Table[ f = TransformedField["Spherical" -> "Cartesian", \[Psi][n, l, m, r, \[Theta], \[Phi]], {r, \[Theta], \[Phi]} -> {x, y, z}]; DensityPlot3D[Abs[f]^2 , {x, -30, 30}, {y, -30, 30}, {z, -30, 30}, ColorFunction -> Hue, ColorFunctionScaling -> True, Boxed -> False, Axes -> False, PlotLabel -> Style[StringForm["Hydrogen atom orbitals\n |\[Psi]\!\(\*SuperscriptBox[\(|\), \\(2\)]\) : n=`` l=`` m=``", n, l, m], Medium, FontFamily -> "DejaVu Serif"], LabelStyle -> {Black, Bold}, RegionFunction -> Function[{x, y, z}, x < 0 || y > 0], PlotLegends -> Automatic], {n, 1, 4}, {l, 0, n - 1}, {m, -l, l}]