Before providing possible approaches:
- Maxwell equations, along with Lorentz force, are the governing equations of classical electromagnetism (treating irregularities as described below);
- Coulomb force directly comes from the combined application of Lorenz force (in steady condition) and Gauss's law for the electric field, when you deal with point charges. It adds nothing to Maxwell equations, and it's a direct consequence.
The differential/local form of Maxwell equations only hold if the fields are "regular enough", being a set of differential equations. If you want/need to use differential form of the equations, I'd try two approaches:
Approach 1. You could "make them work" using generalized functions, like Dirac's delta and steps to treat point charges, surface charges,...
Approach 2. You could 1. partition the domain in a set of regions where the fields are "regular enough"; 2. connect these regions with boundary/jump conditions across the manifolds where "irregularities" arise
Global/integral form of the equations work seamlessly even with functions not regular enough for the differential form to hold, but they usually provide you less detailed information, if you can't exploit symmetry and apply them on particular domains that give you all the information you need. These equation, for a steady domain $V$, read
$$\begin{cases}
\oint_{\partial V} \mathbf{d} \cdot \mathbf{\hat{n}} = Q^{int}_{V} \\
\oint_{\partial S} \mathbf{e} \cdot \mathbf{\hat{t}} + \frac{d}{dt} \int_{S} \mathbf{b} \cdot \mathbf{\hat{n}} = 0 \\
\oint_{\partial V} \mathbf{b} \cdot \mathbf{\hat{n}} = 0 \\
\oint_{\partial S} \mathbf{h} \cdot \mathbf{\hat{t}} - \frac{d}{dt} \int_{S} \mathbf{d} \cdot \mathbf{\hat{n}} = I_{S} \\
\end{cases}$$
being $Q^{int}_V$ the total electric charge in the volume $V$ and $I_S$ the electric current through surface $S$.
Edit. Here you can find a possible approach, using two different models of the thin distribution of electric charge (with no detail known, except for radial symmetry and total charge), one with Dirac's delta and one with uniform charge in a finite-thickness layer close to the surface of the sphere. As you can see, if you use the Dirac's delta model, the electric field is discontinuous at the surface, so it's meaningless to look for a value there (if you don't further regularize the integral of the Dirac delta, finding the average of the values on the sides of the jump); otherwise, if you use the finite approximation, the electric field is continuous and you can find a value of the electric field for each value of $r$.
Edit 2. - Evaluation of the electric field using Gauss' law and generalized function. In this paragraph, I'll show you how to use a differential equation with generalized functions, as the one required to represent a surface distribution in a 3-dimenisonal domain.
Using spherical coordinates, and the assumption of spherical symmetry, here the charge density is represented by the generalized function
$$\rho(r) = \sigma \delta(r-R) \ .$$
Gauss' law $\nabla \cdot \mathbf{e} = \frac{\rho}{\varepsilon}$ in spherical coordinates reads
$$\frac{1}{r^2} \frac{d}{d r} \left( r^2 e_r \right) = \frac{\rho}{\varepsilon} \ ,$$
and it can be recast as
$$\frac{d}{d r} \left( r^2 e_r \right) = r^2 \frac{\rho}{\varepsilon} \ .$$
Direct integration of the last equation reads
$$\begin{aligned}
r^2 e(r) - r_0^2 e(r_0) & = \frac{\sigma}{\varepsilon} \int_{r_0}^{r} x^2 \delta(x - R) d \, x = \\
& = \frac{\sigma}{\varepsilon} R^2 u(r-R) \ ,
\end{aligned}$$
being $u(r)$ the step function
$$u(r) = \begin{cases} 0 & , \quad r < 0 \\ 1 & , \quad r > 0\end{cases}$$
Remark. This is the discontinuous function that appears in the diagrams above, and that you can regularize. But if you use this "Dirac's delta"-approach this function is discontinuous and, as I've already mentioned, it's meaningless to look for a value in $r=R$; if you really want, you can set $\frac{1}{2}$ keeping in mind that it's just the average value of the jump.
Setting $r_0 = 0$, the analytic expression (piece-wise continuous) of the radial component of the electric field reads
$$e(r) = \frac{\sigma}{\varepsilon} \frac{R^2}{r^2} u(r-R) =
\begin{cases} 0 & , \quad r < R \\
\frac{\sigma}{\varepsilon} \frac{R^2}{r^2} &, \quad r > R \ .\end{cases} $$