Landau diamagnetism takes place because of the noncommutativity of
$\mathbf{p}$ and $\mathbf{A}(\mathbf{r})$. In the classical treatment no such noncommutativity exists, thus the magnetic susceptibility is identically zero.
One way to appreciate the dependence of the result on $\hbar$ and at the same time perform the full quantum treatment is to perform the computation in the coherent state basis. In this basis, the Landau Hamiltonian has the form of an isotropic two dimensional harmonic oscillator with the Larmor frequency as a natural frequency, (for some detail, please see the following question .
$H = \hbar \omega_c (a^{\dagger}a + \frac{1}{2})$
Due to the nonvanishing commutation relation between the creation and
anihilation operators, the exponential of the Hamiltonian operator is given by:
$\exp(-\beta H) = \exp(-\frac{1}{2}\beta \hbar \omega_c) \exp(-\beta \hbar \omega_ca^{\dagger}a(1-e^{-\beta \hbar \omega_c}))$
In the coherent space basis, the partition function is given by:
$Z = \omega_c \int d^2\alpha \exp(-\frac{1}{2}\beta \hbar \omega_c) \exp(-\beta \hbar \omega_c\bar{\alpha}\alpha(1-e^{-\beta \hbar \omega_c})) = \omega_c (1-e^{-\beta \hbar \omega_c})^{-1}$
(The multiplicative factor proportional to $\omega_c$ is the Jacobian of the transformation of the phase space volume). (I am being careless in the immaterial constant terms).
This partition function gives the correct magnetic susceptibility. In addition, In the small \hbar limit the partition function becomes a constant, thus giving the Bohr-van Leewen theorem.
...it makes no sense to calculate any property of the particle with the "initial" wave-function, since this is simply the incorrect wave-function for the new well?
The wavefunction can't be "incorrect for the well". Your wavefunction is just an initial condition for time-dependent Schrödinger equation. Here's how it would evolve if you solve the time-dependent equation (I ignore normalization here):
We see that since the actual particle's wave-function is not defined when $x>L$...
The wavefunction is defined for all $x\in\mathbb R$. It's just zero outside of $(0,L)$ because the potential is infinite there.
the second term will be zero
This remains true, however, because of what I said above.
the wave-functino is normalized for $0<x<L$
Actually, again, the wavefunction is normalized period. It's defined for the whole real line, and zero outside the well, so when you normalized using the integral over the well, it's the same as if you integrated over $\mathbb R$. If it were not this way, your "normalization for some domain" would not make any sense, i.e. it'd not be a normalization at all.
Your further calculation looks OK to me.
I think maybe my difficulty with "visualzing" the problem is that I do not fully understand the expression for $c^2_n$, and how this gives the probability for measuring the energy level $n$.
That this gives the probability for measuring energy level $n$ is known as Born rule. You find projection of your actual wave function onto the eigenstate of energy, namely on state $n$. By Born rule, its magnitude squared is the probability of measuring the system to appear in that eigenstate.
The fact that you used only the unchanged original wavefunction for calculations, even though it almost immediately drastically changes as time goes on, is that despite its change in shape, its coefficients $c_n$ actually only change their phase like
$$c_n\propto\exp\left(-\frac i\hbar E_nt\right),$$
but remain the same in magnitude — because the potential is independent of time. Thus you can measure the energy after some time passes, and will still get the same result.
Best Answer
OP is considering the Bohr-Sommerfeld quantization rule
$$ \oint k(x) \mathrm{d}x ~=~2\pi (n + \frac{1}{4}\sum_i\mu_i) , \qquad n\in\mathbb{N}_0,\tag{1} $$ where the sum $\sum_i$ is over turning points $i$ with positions $x_i$ and where $\mu_i\in\mathbb{Z}$ is the metaplectic correction/Maslov index of the $i$th turning point. A turning point comes in different types:
A turning point of generic type has Maslov index $\mu_i=1$. This is seen from the semiclassical connection formulas. The connection formulas with the well-known $\exp(\pm i\frac{\pi}{4})$ phase shift [which induces a $2\times \frac{\pi}{4}=\frac{\pi}{2}$ phase shift between the left and the right mover, and corresponds to $\mu_i=1$] are derived under the assumption that there exists an overlapping region where the following two conditions are both satisfied:
The quasi-classical condition: $|\lambda^{\prime}(x)| \ll 1$. This typically holds away from the turning point.
A linearization of the potential $V(x)$ is valid. This typically holds only in a small neighborhood around the turning point.
An infinite square well potential (say located at $x=x_i$) typically satisfies condition 1 for $x\neq x_i$ but fails condition 2. The connection formula is then replaced by a boundary condition (BC) $$ \psi(x_i)~=~0,\tag{2}$$ cf. e.g. this Phys.SE post. The BC (2) implies a $\pi$ phase shift between the left and the right mover, which corresponds to a Maslov index $\mu_i=2$.
For more details, see also e.g. the Einstein-Brillouin-Keller (EBK) quantization rule and this Phys.SE post. See also this Phys.SE post for references.