For a better description of what follows, see Sec. 4.2.3.3 of A Guide to Monte Carlo Simulations in Statistical Physics by the giant Kurt Binder.
The punchline is that it depends on whether you sample two phases or not. Usually, the finite susceptibility is the one measured by sampling only one phase. When sampling two phases, such as your case with $h\to 0$, the measured susceptibility would diverge. In the thermodynamics limit, it is a matter of which limit you take first. If you take first $h\to 0$ and then $L\to \infty$ then you are not putting the system in a given phase therefore, the system can sample two phases. However, if you set $L\to\infty$ and then $h\to 0$, you're blocking your system into one given phase.
To understand that, I'll define a probability distribution of single spin $s$ in a mixture:
$$P_L(s)=A\left(e^{-\displaystyle{L^d}\left[\dfrac{(s-m)^2}{2\chi_{\rm pure} T}-\dfrac{hs}{T}\right]}+e^{-\displaystyle{L^d}\left[\dfrac{(s+m)^2}{2\chi_{\rm pure} T}-\dfrac{hs}{T}\right]}\right),$$
away from the critical point, this should be a good description. What it tells you is that, for an ensemble of system. You have a bimodal distribution of spin $s$ around $\pm m$ (the mean field solution of ising for example) with a single phase susceptibility $\chi_{\rm pure}$ (because you are just adding two "independent" phases). On top of that, you add a magnetic field $h$ making things aligned.
From this you can compute the susceptibility of the mixed phase (which is equivalently given by the variance of $s$ from $P_L$ from the fluctuation dissipation theorem):
$$\chi_{\rm mixed}=\dfrac{\partial\langle s\rangle}{\partial h}=\beta L^d\left(\langle s\rangle^2-\langle s\rangle^2\right)=\chi_{\rm pure}+L^d\dfrac{\beta m^2}{\cosh^2\left(\beta hm L^d\right)}.$$
Now, you find the interesting limits:
$$\lim_{L\to \infty}\lim_{h\to 0}\chi_{\rm mixed}=\chi_{\rm pure}+L^d\beta m^2,$$
which diverges due to the bimodality of the distribution, while:
$$\lim_{h\to 0}\lim_{L\to \infty}\chi_{\rm mixed}=\chi_{\rm pure}$$
does not diverge and converge to the value of the susceptibility found for the pure phase. As explained above, this is understood from the fact that $\lim_{h\to 0}\lim_{L\to \infty}$ selects a single phase by symmetry breaking while for $\lim_{L\to \infty}\lim_{h\to 0}$ both phases are equally sampled. This last case corresponds to your definition (because you take into account the $\uparrow$ and $\downarrow$ states) and hence the divergence you observe. With your computations, a way to obtain a single phase would be to restrict yourself at $h\to 0^\pm$.
Concerning the scaling, the above works far from the critical point due to the Gaussian assumption. Closer to the critical point, the renormalization group provides an other answer. It works for both the mixed state or the single phase susceptibility.
The (singular part of the) susceptibility must be a homogeneous function of $h$, $t$ and also the system size: $L$ (which sets another relevant length):
$$\chi(t, h, L)\propto L^{-\gamma/\nu}\mathcal F(tL^{1/\nu}, hL^{(\beta+\gamma)/\nu}),$$
with $\mathcal F$ a non-universal function.
The exponents are chosen such that the correct exponents are recovered in the thermodynamics limit. For example at $h=0$ and for large system size, in the critical region, we must obtain $\chi(t, h=0, L\to\infty)\propto t^{-\gamma}$. Indeed, we obtain this scaling (the only possible choice to eliminate the $L$ dependence being $\mathcal F(u, 0, \infty)\propto u^{-\gamma}$):
$$\chi(t, 0, L\to\infty)\propto L^{-\gamma/\nu}(tL^{1/\nu})^{-\gamma}\propto t^{-\gamma}.$$
It also tells you how the susceptibility grows at the critical point:
$$\chi(t\to 0, 0, L)\propto L^{-\gamma/\nu},$$
which is infinite at infinite $L$. For the mixed susceptibility, you will also grow infinitely with this scaling below the critical point for increasing system size.
