Why is the potential not infinite? One way to calculate potential (using infinity as our reference point) is to sum all the contributions of charges that are around.
Let's say I want to calculate the potential at some point on charged surface.
At that point, there is some charge (can be infinitesimal) and that charge should contribute something divided by zero (since the distance is zero) to the potential at that point. Using that logic, every point on charge distribution should have infinite potential.
What is flawed with this argument?
 A: No, it would not necessarily mean infinite. This is a classic mathematical misunderstanding regarding limits.
If you have a fraction,
$$\frac ab,$$
and you let the numerator tend to zero, $a\to 0$, then the fraction might tend towards zero:
$$\frac ab\to\frac 0b=0\; \text{ for } \;a\to 0.$$
If you instead let the denominator tend to zero, $b\to 0$, then the fraction might tend towards infinity:
$$\frac ab\to\infty\; \text{ for } \;b\to 0$$
But what if both happen simultaneously? Which one will then win? Is the numerator or the denominator winning? Will the fraction as a whole tend towards zero or infinity, or something in between (converging)?
This is the situation you have. And you have entirely dismissed the numerator which is infinitesimal and are then claiming that the whole fraction is only based on the denominator tending towards zero. This is a mathematical misunderstanding. Firstly we can't directly know the answer, and secondly the answer depends on "how much" or "how fast" the number in the numerator and in the denominator tend towards their limits.
A: Physics equations are models (or abstract descriptions) of observed behavior - and as such they make certain assumptions, and so their domain of applicability covers those scenarios where those assumptions are valid.
Mathematically, when the distance is zero, you have a zero in the denominator, and the value at that point is, technically, not infinite, it's undefined - there's, as mathematicians would say, a singularity there, due to discontinuity, and the model doesn't apply. For any other point arbitrarily close to zero, the value can get large (and tends to infinity as you approach zero), but it's finite everywhere.
But that's not the end of the story - if you want to treat some small charged volume as a dimensionless point, then when you get close enough, the charge stops being point-like (so the point-charge assumption breaks down), and on an even smaller scale, quantum mechanical effects become important.
Now, back in the macroscopic realm, if you're treating a charge distribution as a continuous charge density, note that an infinitesimal volume (or surface) element does not have zero volume (or zero surface) - the notion is just that it's "very small". Charge density is also a model, a way to treat things as a continuum, and ignore the underlying discontinuous nature (as it doesn't really come into play at the scales you're working at). Conceptually, charge density at a point doesn't quite mean that there's charge at that (mathematical) point; rather, it means that in the immediate (infinitesimal) surroundings of that point, you can take the charge density to be constant, and recover the charge for a small volume by multiplying the volume with the density.
Pretty much in the same way that a mass density at a mathematical point doesn't mean there's mass there (I'm switching to mass density because IMO it's more familiar and potentially easier to grasp). E.g., consider a gas cloud that varies in density throughout its extent. If you pick an arbitrary point in a gas cloud, you are likely to pick one that's between the molecules (i.e. there's nothing there, the point doesn't actually contain mass). But as you expand a volume around that point, you engulf more and more molecules (a huge number of them in quite a small volume), and if the volume is small enough, the density of the gas in it can be treated as roughly uniform. That's what density models (describes abstractly) - at macroscopic scales.
A: 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.
A: Point charges are just an abstraction, a limit model, like surface charge distribution in a infinitely thin layer.
You can think at a point charge as the limit for $r \rightarrow 0$ of a continuous distribution in a volume. As you look inside the small charges, the charge enclosed in the volume decreases and goes to zero when you're looking to a point.
The actual field depends on the distribution of the charges inside the small volume. If we assume a uniform charge distribution inside the spherical volume, we end up with an electric field
$\mathbf{e}(\mathbf{r}) = \left\{ \begin{matrix} \frac{q}{4\pi\varepsilon}\frac{\mathbf{r}}{R^3}\quad r \le R \\ \frac{q}{4\pi\varepsilon} \frac{\mathbf{r}}{r^3} \quad r \ge R \end{matrix}\right.$
and potential $\phi (\mathbf{r}) $ s.t. $\nabla \phi(\mathbf{r}) = \mathbf{e}(\mathbf{r})$
$\phi(\mathbf{r}) = \left\{ \begin{matrix} \frac{q}{8\pi\varepsilon}\frac{r^2}{R^3}-\frac{3}{8}\frac{q}{\pi\varepsilon}\frac{1}{R}\quad r \le R \\ -\frac{q}{4\pi\varepsilon} \frac{1}{r} \qquad \quad \quad r \ge R \end{matrix}\right.$
where the constant $\frac{3}{8}\frac{q}{\pi\varepsilon}\frac{1}{R}$ (if I made no mistake) appears to get a continuous potential field going to zero at infinity (this condition at infinity sets the value of the irrelevant additive constant for our choice of the potential).
You can use other distributions, with the given total charge. The results change only locally inside the charge itself. But if you use a model of point charges your main interest is not the distribution inside the small volume, and thus every distribution inside the small mass with the same total charge should be equivalent for your model and your goal.
The same applies to surface distributions, and other singular charge distributions. You can be rigourous mollifying the differential equations to introduce point charges to regularize delta-functions.
If you wish to start with the Maxwell's equations as the fundamental equations of electromagnetism, a point charge of intensity $q$, located in $\mathbf{r}_0$ is represented by a charge density
$\rho(\mathbf{r}) = q \delta( \mathbf{r} - \mathbf{r}_0)$
A: We are treating the charge as a volume charge density. Now you say that at any point why isn't the field infinite.
It isn't so because exactly at that point the charge is zero. Why? Because a point has no volume and hence charge at that point is zero. Hence
$Q = 0$
when
$\Delta v = 0$
$\Rightarrow$
$\overrightarrow{\mathrm{E}} \neq \infty$
