Quantization of electrons' angular momentum in atoms and molecules It is known that the Schrödinger's equation of the electron's wave function in atoms can be solved analitically only when a single electron is present (the "hydrogenlike atom"). In that case, the angular part of the equation leads to the solution with spherical harmonics, whence the quantization rules and numbers for angular momentum.
When speaking about atoms and molecules with more than one electron, it is widely assumed that the angular momentum is quantized in a similar manner. This is a good approximation based on experimental observations, or it can be proved mathematically that the azimuthal and magnetic quantum numbers can be used notwithstanding the presence of several electrons and in molecules of several nuclei (in these cases the Schrödinger's equation does not contain a potential function that is purely radial)?
According to an answer to a previous question, it seems that it is an observational evidence rather than a mathematical assertion (https://physics.stackexchange.com/q/22815).
 A: This kind of depends on exactly what it is you're talking about.


*

*In multi-electron atoms, the total angular momentum $\mathbf J = \mathbf L+\mathbf S$, which includes the orbital and spin angular momenta of all the electrons in the atom, is rigorously conserved. This is the $J$ quantum number of the atom and it is always a "good quantum number" (in the sense that it represents a conserved quantity).

*If you ignore spin, and spin-orbit coupling, then even in a multi-electron atom the orbital angular momentum is rigorously conserved, which seems to be one of your primary worries when you say things like

in that case the potential is not central, because of the repulsion between the different electrons.

Indeed, for multi-electron atoms the potential is no longer central, but it is still symmetric under global rotations: if you move a single electron without moving the others then the hamiltonian obviously changes (so the angular momentum of each individual electron isn't conserved), but if you move them all under the same transformation then their relative distances are untouched, so the interaction remains invariant. This global symmetry under rotations means that the generator of that symmetry (i.e. the total orbital angular momentum operator, $\mathbf L = \sum_j \mathbf L_j$) is rigorously conserved, and there is a shared eigenbasis of the electronic hamiltonian and this angular momentum.

*In multi-electron atoms, the total orbital and spin angular momenta (i.e. $\mathbf L$ and $\mathbf S$ separately, each of which includes contributions from all the atom's electrons) are rigorously conserved in the absence of strong spin-orbit coupling. This is the framework is where term symbols of the form $^{2S+1}L_J$ with well-defined orbital and spin angular momenta (like, say, ${}^2\mathrm{P}_{3/2}$) come from.
As the strength of the spin-orbit coupling increases, which generally happens as the atomic number gets bigger, you will generally start getting a small amount of mixing between terms of different character (so, say, a ${}^2\mathrm{P}_{1/2}$ state might also include a 0.1% amplitude of ${}^4\mathrm{D}_{1/2}$ character). In more rigorous language, this means that the hamiltonian eigenfunctions no longer have 100% overlap with their dominant contributing state in a definite-$\mathbf L$ basis, or that the $\mathbf L$ eigenfunctions are no longer exact eigenfunctions of $H$.
(Nevertheless, though, the fact that this involves an approximation on the side of which terms contribute to the hamiltonian should not be taken as a way to immediately dismiss angular-momentum conservation. In particular, everything in physics is approximate, and you always need to look at what other approximations are involved in whatever you're doing. In this specific case, if you wanted to detect that ${}^4\mathrm{D}_{1/2}$ contribution, then you'd look for a transition line that would otherwise be (dipole-)forbidden in the ${}^2\mathrm{P}_{1/2}$ state ─ like, say, to an $\rm F$ state ─ but then you have to seriously consider whether an observation of that transition is evidence of that ${}^4\mathrm{D}_{1/2}$ contribution or whether it is instead caused by quadrupole or higher-order selection rules.)
As the atomic number gets really big, and you get into the large-atom regime, then $\mathbf L$ stops being a "good quantum number" (i.e. the spin-orbit coupling becomes so big that $\mathbf L$ eigenfunctions are no longer good approximations to the hamiltonian eigenstates), and you pass through some intermediate coupling schemes (which again have mostly-well-defined angular momenta coming from sundry combinations of spin and orbital angular momenta from different shells) before settling into jj coupling at the high-$Z$ end of the scale, where you have states with well defined $\mathbf j$ per shell (with trace contributions from other $j$'s).

*In multi-electron atoms, it is also common to talk about the angular-momentum characteristics of individual electrons, and the orbitals they occupy to form the system's multi-electron state. These orbitals should only ever be considered as a convenient basis within which to formulate the true multi-electron state, and indeed single-electron orbitals have extremely limited physical reality in a multi-electron setting, as I've explained elsewhere.

*For molecules, and specifically for linear molecules like diatomics, we often talk about the angular momentum of the electrons. Generally, the total angular momentum of the electrons is never conserved in such molecules, since the hamiltonian is not spherically symmetric, but within the Born-Oppenheimer approximation we do have one axis of symmetry and therefore the conservation of one component of the angular momentum, generally taken on the $z$ axis. This is where molecular term symbol notation comes from, with one letter denoting the orbital angular momentum (now in Greek, to emphasize that it's only $L_z$ that's conserved; as in atoms, with lower-case terms denoting single-electron orbitals and upper-case terms denoting multi-electron combinations).
This is typically a rigorously-conserved quantity, unless you're looking at a system where nonadiabatic beyond-Born-Oppenheimer couplings between nuclear and electronic motion are non-negligible. (Hint: if you have to ask, then those couplings are negligible.)

*For molecules, we also care about the total angular momentum of the overall motion of the system, i.e. of the nuclear positions. This is a rigorously conserved quantity (so long as you're in the gas phase, obviously), giving a guaranteed "good quantum number" for rotational molecular spectroscopy; any experimental limits on this are purely on the instrumentation side - i.e. this is always conserved, but if the system is too large then that might be hard to verify explicitly.
I'm unsure about what's the current record in terms of the molecular size and mass at which rotational spectroscopic verification of the conservation of this angular momentum is still feasible.
On the other hand, for small molecules, this quantity can be pushed to ridiculous extremes in terms of the value of $J$, which can be pushed all the way up to the hundreds and even a thousand using optical centrifuges (examples: classic, modern).
A: As an initial remark, the Schrodinger equation can be solved exactly for a variety of potential, not just “hydrogen-like” atoms - the cases of the harmonic potential, the Morse potential, or Poschl-Teller potential immediately come to mind but there are multiple other ones as well that don’t have entries on Wikipedia.
Quantization of angular momentum follows exactly from the commutation relation of the angular momentum operators $\hat L_x,\hat L_y$ and $\hat L_z$. These commutation relations do not depend on the number of particles in the system.
The simplest example would be two spin-1/2 particles, made famous by various versions of Bell’s theorem.  The total spin (or angular momentum) of this 2-particle system remains quantized and it can only take the values $0$ or $1$ depending on how the state is prepared.
Another example is the nuclear $su(3)$ model - not so used in chemistry but still quite useful for multi-particle nuclei.  The 3-dimensional harmonic oscillator can be solved using the $su(3)$ Lie algebra, and the angular momentum operators are in this algebra.  It is also “easy” to construct multiparticle $su(3)$ state using Lie algebraic techniques (just an extension of Clebsch-Gordan technology for angular momentum).  For these multiparticle states, angular momenta of the individual constituents are combined in the usual way and remain quantized.
More generally, when the potential is not central, angular momentum will not be conserved for individual single particle states and so states will not necessarily be eigenstate of angular momentum.  An example of this is the Nilsson model for deformed nuclei.  This does not mean that basis states with good angular momentum quantum number cannot be used to start the calculations, just that the final states will be a linear combination of those single particle basis states and that $\Delta \ell\ne 0$.  
Of course since the total Hamiltonian must be a rotation scalar, the final eigenstates of $H$ can be chosen to have good angular momentum.  Angular momentum remains quantized since the total angular momentum operator still satisfy the usual commutation relations.
A: To sum up:


*

*the total angular momentum $\vec{J}$ is always conserved in isolated atoms and molecules, since the potential functions do not depend on a particular direction (they are rotational invariant for a global rotation of the system around the origin of the frame of reference, tipically the center of mass for multi-nuclear molecules) and thus the wave function of the system is an eigenfunction of the angular momentum squared magnitude operator $\hat{J}^2=\hat{J}^2_x+\hat{J}^2_y+\hat{J}^2_z$ ($\hat{J}_i$ is the generator of infinitesimal global rotations around the axis $i$).

*since it is possible to construct mathematically a angular momentum operator for a single direction $z$ that commutes with the operator $\hat{J}^2$, namely $\hat{J}_z$, also the angular momentum projection $J_z$ is conserved (see e.g. Landau, Lifshitz,  Course of theoretical physics, III vol., § 27).

*from the conservation of $J_z$ and the mathematical expression of the operator $\hat{J}_z$, it is possible to define a total angular momentum projection quantum number $M_J$, so that $J_z = M_J \hbar$. The number $M_J$ must be integer or semi-integer (in order to have a single-valued probability density). Also the squared magnitude of the total angular momentum can be derived as $\vec{J}^2=J(J+1)\hbar^2$, where $J \geq |M_J|$ is the total angular momentum quantum number, also integer or semi-integer.

*the orbital angular momentum of all the electrons $\vec{L}$ and the spin angular momentum of all the electrons $\vec{S}$, whose sum is the total angular momentum $\vec{J}=\vec{L}+\vec{S}$ (disregarding the nuclear angular momentum $\vec{I}$) are approximately conserved separately, i.e. they are conserved separately if the influence of the spin is disregarded in the solution of the Schrödinger equation of the electrons' wave function (if $\vec{S}$ is fixed, the conservation of $\vec{L}$ follows from the conservation of $\vec{J}$). In that approximation, it is possible to define the azimuthal quantum number of the orbital part of the electrons' angular momentum $L$, the magnetic quantum number of the orbital part of the electrons' angular momentum $M_L$, the spin quantum number $S$ and the spin projection quantum number $M_S$, that are good quantum numbers.

*it is not possible to separate a conserved angular momentum of the orbit of the single electrons in atoms and molecules, so $\ell$, $m_\ell$ are not good quantum numbers. Also the spin projection quantum number $m_s$ is not conserved separately for any electron (of course, the spin quantum number is always $s=1/2$ because that is an intrinsic property of electrons). However, in single atoms the approximation of electrons staying in separate orbital with definite $\ell$, $m_\ell$, $m_s$ is a good picture (at least for light atoms).

*the same is a fortiori true for molecules, where additionally it is more difficult to separate a $L$ term, since the orbital angular momentum must be calculated with reference to the center of mass, not to a single nucleus. Therefore, the orbitals with a specific magnetic quantum number are defined only for simple molecules (e.g. diatomic molecules with an axial symmetry).

