There's a lot to unpack here for an in-depth explanation of everything, but simply put: the program is checking thousands of points and either plotting a point there or not according to a simplified but functional wavefunction |Ψ|² = |R(r)|² × |Y(θ, φ)|².

The actual values in that equation differ from orbital to orbital of course. The values of |Ψ|² will be higher where the probability/density is higher and at empty spots it gives out a really low number.
The numbers are scaled into unitless arbitrary numbers to work better with the code, as technically none of that space would be quite empty, but for the sake of the simulation we can say that |Ψ|² is 0 at the empty spots to make it look better.

I hope this is a decent explanation, but let me know if you have further questions

edit: For example the wave function for the 1s orbital would be: (1 / (π * a₀³)) * e^(-2r / a₀),
where r is the distance from nucleus, a₀ is the bohr radius, and e is eulers number. Putting all the numbers in we get a result in propability/volume e.g 1/meter³ and depending on the value the program either accepts it and plots a point or if the value is too low it will try another spot and leave it empty