So to make that more explicit. Take the case of the uniformly charged ball:
sergiokapone said:
Let's say, for example, we can take the problem of finding the field of a uniformly charged ball. Here it is obvious that integral equations and symmetry are enough. But if you solve differential equations, you would need to know how to glue the solutions together and define constants.
... or actually ... let's spice it up ... uniformly charged spherical
shell so we actually have some delta distributions in the PDE. The field is then described by the differential equation$$
\nabla \cdot \vec E = \frac{Q}{4\varepsilon_0\pi R^2} \delta(r- R).
$$ Here ##R## is the radius of the sphere. The spherical symmetry of the problem results in ##\vec E = E(r) \vec e_r## and so $$
\nabla \cdot \vec E = \frac{1}{r^2}\frac{d(r^2 E(r))}{dr} = E'(r) + \frac{2E(r)}{r} = \frac{Q}{4\varepsilon_0\pi R^2} \delta(r- R).
$$ This is now an ODE that we need to solve. For brevity of notation, let's introduce ##C = Q/4\varepsilon_0\pi R^2##. Since the inhomogeneity on the RHS contains a delta function at ##r = R## and ##\delta(r-R) = \theta'(r-R)##, we make the ansatz ##E(r) = f(r) \theta(r-R) + g(r) \theta(R-r)##. Clearly ##f(r)## describes the behaviour of ##E(r)## for ##r > R## and ##g(r)## for ##r < R##. Both ##f## and ##g## are assumed to be smooth functions although it really does not matter what ##g(r)## is for ##r > R## etc. The derivative ##E'(r)## is now given by
\begin{align*}
E'(r) &= f'(r)\theta(r-R) + f(r) \delta(r-R) + g'(r) \theta(R-r) - g(r)\delta(r-R) \\
&= f'(r) \theta(r-R) + \delta(r-R)[f(R) - g(R)] + g'(r) \theta(R-r).
\end{align*} Insertion into the differential equation leads to $$
[f'(r) + 2f(r)/r] \theta(r-R) + \delta(r-R)[f(R) - g(R)] + [g'(r) + 2g(r)/r] \theta(R-r) = C \delta(r-R)
$$Three things are necessary to satisfy this:
\begin{align*}
f(R) - g(R) &= C \\
f'(r) + 2f(r)/r &= 0 \quad (r > R) \\
g'(r) + 2g(r)/r &= 0 \quad (r < R).
\end{align*}
For ##f## and ##g## we have the same form of differential equation, resulting in the same functional form of the general solution:
$$
f(r) = \frac{B_f}{r^2}, \quad g(r) = \frac{B_g}{r^2}.
$$
For ##g## we must require that ##B_g = 0## in order to avoid the resulting field ##\vec E(r)## having a delta divergence at ##r = 0##*. We now have ##f(r) = B_f/r^2## and ##g(r) = 0##. This implies that $$
f(R) - g(R) = f(R) = \frac{B_f}{R^2} = C = \frac{Q}{4\varepsilon_0\pi R^2}$$ and therefore $$
B_f = \frac{Q}{4\varepsilon_0 \pi} \quad \Longrightarrow \quad \vec E = \frac{Q}{4\varepsilon_0 \pi r^2} \theta(r-R) \vec e_r$$
The above essentially derives Newton's shell theorem. The field for a uniformly charged ball can be obtained by integrating over all infinitesimally thick spherical shells out to the ball radius ##R_1##, noting that the charge on each is ##Q = 4\pi \rho R^2 dR##, where ##\rho## is the constant charge density: $$
E_{\rm ball}(r) = \frac{\rho}{\varepsilon_0 r^2} \int_0^{R_1} \theta(r-R)R^2 dR =
\begin{cases}
\frac{\rho R_1^3}{3\varepsilon_0 r^2} & \qquad (r > R_1) \\
\frac{\rho r}{3\varepsilon_0} & \qquad (r < R_1)
\end{cases}
$$
Note that ##\rho = 3 Q_{\rm tot}/4\pi R_1^3## and so this may also be written$$
E_{\rm ball}(r) =
\begin{cases}
\frac{Q_{\rm tot}}{4\pi \varepsilon_0 r^2} & \qquad (r > R_1) \\
\frac{Q_{\rm tot} r}{4\pi \varepsilon_0 R_1^3} & \qquad (r < R_1)
\end{cases}
$$ which may be more familiar.
Further notice that the integral form of the field equations are nowhere in sight. In fact, you can easily derive the integral form from the differential form (and vice versa). In the differential form to integral form direction:
Assume that ##\nabla \cdot \vec E = \rho/\varepsilon_0##. Then for any volume ##V## with boundary surface ##S##, the flux of the electric field through ##S## is given by $$
\Phi = \oint_S \vec E \cdot d\vec S$$ By the divergence theorem, we therefore find that$$
\Phi = \int_V \nabla\cdot \vec E \, dV = \frac{1}{\varepsilon_0} \int_V \rho\, dV \equiv \frac{Q_{\rm enc}}{\varepsilon_0}$$ where ##Q_{\rm enc}## by definition is the total charge in the volume ##V##.
In true textbook fashion, I will leave the other direction for the interested reader.
* Note that ##r = 0## is not really a boundary of the physical region in which we are solving the PDE! It is common to argue that ##r = 0## is a boundary and that hand-wavingly state that the solution must be finite in that regime. That is not a very good argument though. The better argument comes from the ##1/r## field having a delta-divergence at ##r = 0## meaning that it does not actually solve the ##\nabla \cdot \vec E## at ##r = 0##.
Edit: Fixed an error in the solution. (I am too used to solving for the potential rather than field strength apparently…)