You can use an Emission shader for the Volume socket of your Material output node, and give it a volumetric texture.
It's pretty easy for a 1s orbital, (n,l,m) = (1,0,0). (spherical symmetry, monotonic decay: it's not actually exponential in my example below):

Create a sphere, and assign a material to it like the one in the picture.
The first three nodes to the left (Geometry, Object Info and Vector Math>Subtract) will output a vector (x-x₀, y-y₀, z-z₀) for each x,y,z point in the sphere's volume, where x₀,y₀,z₀ is the pivot point of the sphere. The following node Vector Mapping can be used to scale the texture. These are the cartesian coordinates that will be used to define a scalar function, like a wavefunction.
The last four nodes to the right ("Math>Absolute, Color Ramp, Emission, Material output) will take the output of your scalar function and assign it to the volumetric emission, colored according to a color map (I've taken the absolute value, you might want to take the 2nd power, or to have different colors for the positive and negative values).
In the middle, you will need your actual scalar function. It can be any node, or group of nodes, that takes a Vector (blue socket) as input, and outputs a scalar (grey socket). In my example, I've used a Gradient texture node. Unfortunately, complex functions are not built-in in blender. You can, however, create a Script node that implements the wavefunctions (you will probably need numpy.polynomial.legendre
)*
* I've never done this before, so I can't guarantee it's possible to just make python do the math. If it's possible, I think this wiki page is a good starting point.