To answer this question, you need some model of how the non-pointlike part of the charge distribution looks. In our three-dimensional universe, nuclear matter has approximately constant mass density and charge density.
The potential energy $V$ of a charge $e$ moving in an attractive potential with this shape is
$$
V_3(r) = -e\frac{q_\text{enclosed}}{r} \sim
\cases{
- {e\rho R^3}{r^{-1}} & outside, $r > R$
\\
-V_\text{center} + e\rho r^2 & inside, $r < R$
}
$$
Here $\rho$ isn’t quite the charge density, because there is some geometrical factor which will be different from $4\pi/3$ in $d\neq 3$ dimensions.
There is also some magic that makes the potential surface smooth if you are careful about doing the integral from infinity down to the interior of the sphere, and may introduce some factor of two or one-third to the harmonic oscillator term.
But that’s beside the point: inside a uniform charged sphere, the Coulomb potential turns into a harmonic oscillator potential.
If the appropriate four-dimensional potential is
$$V_4(r) = -e\frac{q_\text{enclosed}}{r^2},$$
as suggested in the linked answer, and the hypervolume of the uniform-density charge region goes like $R^4$, you should still recover the harmonic oscillator potential inside:
$$
V_4(r) \propto \cases{
r^{-2} & outside, $r > R$
\\
r^{+2} & inside, $r < R$
}
$$
It’s my recollection that the harmonic oscillator potential separates, $r^2 = x^2 + y^2 + z^2 + w^2$, and so you’ll have a spectrum of partly-degenerate bound states with energies $\hbar\omega(n_x + n_y + n_z + n_w + \frac42) - V_\text{center}$.
Those energies should be relatively un-perturbed so long as the expectation value for the wavefunction outside of the charged region is small.
It might be amusing to examine the correspondence between these harmonic-oscillator states and the stability criteria described in the linked answer and its references.
You ask in a comment about multi-electron states. Recall that the previously-linked answer defines the dimensionless parameter
$$ Z = \frac{2me_{d=4}^2}{\hbar^2}$$
with the results
- $Z \leq 0$: no bound states, i.e. ionized atom
- $0 < Z \leq 1$: the spectrum may or may not be bounded from below
- $1 < Z$: energy spectrum is definitely unbounded from below
Those definitions are for the hydrogen atom, where the nucleus and the electron have the same charge magnitude $|e|$. For this answer we have written $eq_\text{enclosed}$ instead of $e^2$. For a pointlike nucleus with charge $N_p e$, and $N_e$ electrons in pointlike bound states, you might say that the $(N_e+1)$-th electron sees
$$
q_\text{enclosed} = (N_p - N_e)e
\\
Z_\text{effective} = \frac{2me^2}{\hbar^2}(N_p-N_e)
$$
That is, if the attraction is “too strong” $Z>1$ to allow a ground state in the $r^{-2}$ part of the potential, but the details at $r=0$ conspire so that there is a ground state somewhere, then additional electrons might see a “shielded charge,” a correspondingly weaker attraction, and have a lowest-energy state in the $r^{-2}$ part of the potential. The extreme multi-electron state with $N_p - N_e \leq 0$ likewise corresponds to $Z \leq 0$: it makes sense that a neutral atom or a negative ion doesn’t have additional hydrogen-like bound states.
Of course, if your intro-quantum book goes as far as helium, you already know that the “shielded charge” model is qualitative and hand-wavy at best. Furthermore this entire discussion is for one particular choice of charge distribution near the origin.
If you have an idea for a non-uniform charge distribution of finite size, you’ll have to do your own integrals about it.
