Short Answer. In classical electrodynamics, point charges are an artificial construct built from the more general charge density. In a charge density, the potential or field does not have an infinite or undefined value, so no problems arise.
Long Answer. In electromagnetics we describe quantities of charge and current in densities $\rho$ and $\mathbf J$ respectively. It is known that when we calculate e.g. the potential of a point charge, the value at the point charge's position is undefined, because the separation goes to zero. One may reasonably ask, then: If every charge represents a discontinuity in $\Phi$, how can a volume charge density have a well-defined value anywhere?
In fact, with classical electrodynamics, we can take charge density to be fundamental, while a point charge is an approximation for an infinitesimal volume containing a fixed amount of charge. The divergence of potential at a point charge might be accepted as a failure of the classical theory of electrodynamics, even if such a model is useful (even considered fundamental) in constructing the theory.
All this being said, let's resolve the paradox. Assume that electric charge is a continuous quantity (accepting the fact that, at small scales, this assumption breaks down; we will return to this later). Coulomb's law may be formulated from the fact that two electrically charged conducting spheres, with total charge (integrated about the surface) of $q_1$ and $q_2$, with centers at $\mathbf r_1$ and $\mathbf r_2$ have an electrostatic force $\mathbf F$ between them, when separated by a distance much larger than either sphere's radius, and the magnitude of this force is:
$$
|\mathbf F| = \frac{1}{4\pi\epsilon_0} \frac{q_1 q_2}{|\mathbf r_2-\mathbf r_1|^2}.
$$
This is determined from experiment (Coulomb, Maxwell, etc). Under the hypothesis that each infinitesimal charge density attracts or repels every other infinitesimal charge density according to an inverse-square law, we find agreement with experiments involving spheres whose radii are comparable to the separation distance (Poisson, 1813).
From this we may build up the theory of electrostatics, and we find that the potential due to a charged sphere in isolation is the same regardless of the sphere's radius; it is a function of the sphere center only, and thus we may treat problems in electrostatics without resorting to finite spheres, by only performing calculations on the sphere centers, and using the assumption that the separation distance is much larger than the sphere radii. But as we attempt to calculate forces when the separation is sufficiently small, this assumption breaks down, and the expression above is no longer applicable.
In this framework, which follows historical lines, the concept of a point charge is derived from the continuous quantity of electric charge. Thus no paradox arises when attempting to extend results with point charges to results with charge densities. The electric field and electric potential within a charge density are not subject to the divergences we find with point charges. It is only by introducing point charges that such problems arise.
Of course, electric charge is in fact quantized, and charge carriers are, to a first approximation, point charges. Although classical electrodynamics can predict the behavior of such point charges as they interact with one-another, it is still subject to restrictions of finite separation; the theory says nothing (and can predict nothing) about what happens when point charge separation approaches zero. Paradoxes of surface or volume charge densities must now be managed by either (a) adapting to a more accurate (quantum) theory, or (b) staying within the framework of classical electrodynamics and treating extremely large numbers of point charges. Let's take the latter approach.
Modeling Discrete Distributions as Continuous. The fundamental idea is that large numbers of point charges behave like a continuous charge distribution as long as we don't shrink the volume too much. To make this idea exact, we will use Gauss's law.
Gauss's law states that the electric flux through a surface is equal to the enclosed charge within that volume, scaled by a constant. This only uses results which were derived for observed charge densities, which appear continuous in most experiments. Thus we have
$$
\oint_{\partial V} \mathbf E \cdot \hat n\ da = \int_V \frac{\rho}{\epsilon_0}\ dV
$$
Where $\partial V$ indicates the surface bounding the volume $V$. The term $\rho$ is the charge density, and we do not assume its form; i.e. it could be continuous or discrete, in either case the above holds. From this equation, we can derive Poisson's equation, which is the differential equation that determines the electrostatic field at every point within the volume. Our goal is to show that for a true charge distribution $\rho$ which is made up of point charges, there exists an equivalent charge distribution $\rho'$ which is continuous and gives the same results, so long as the volume $V$ is much larger than the distance between point charges.
Suppose we have $\rho$ as the sum of a large number of point charges $\{q_i\}$, and by some process we find an equivalent continuous distribution $\rho'$ such that
$$
\sum_V \frac{q_i}{\epsilon_0} = \int_V \frac{\rho'}{\epsilon_0}\ dV
$$
For a given volume $V$. When the volume $V$ is perturbed, the sum on the left changes by including additional charges, or excluding charges originally present. The density $\rho'$ should be chosen so that an equal change occurs when integrating the continuous distribution, as long as the change in volume includes or excludes a large number of the point charges. This way, whatever the volume $V$ becomes, the expression above holds, so that the electric flux through the surface (which depends only on the total enclosed charge) is the same whether we use the point charges $\{q_i\}$ or the continuous distribution $\rho'$.
Thus far we have not shown that the distributions are equivalent (they are not), but we have shown that the integrated behavior of $\mathbf E$ is nearly the same regardless of which distribution we choose. Now we will apply the divergence theorem to the left-hand side of Gauss's law in integral form, and we get
$$
\int_V (\nabla \cdot \mathbf E)\ dV = \int_V \frac{\rho}{\epsilon_0}\ dV,
$$
and we observe that if this holds for every volume $V$, then we can generally assign
$$
\nabla \cdot \mathbf E = \frac{\rho}{\epsilon_0}.
$$
This is Poisson's equation. When $\rho$ is composed of point charges, it is in fact true that the equation above holds for any volume $V$, so Poisson's equation can be used with the true distribution $\rho$.
It is not true that it holds when we use the approximate continuous distribution $\rho'$, because for very small volumes or even very small changes in volume, our approximation of a continuous distribution breaks down. However, we can still use Poisson's equation under the conditions where our approximation is accurate! Thus we can assign an exact value of the electric field $\mathbf E$ within the charge distribution $\rho'$, despite the fact that the true electric field in that volume will not be what our calculations indicate. The integrated behavior, however, will be sufficiently precise, and this is what we intended to show.
