Can $E=\frac{q}{4\pi\epsilon_0 r^2}$ be directly derived from differential form of Maxwell equations? The electric field of a point charge $q$ is well known to be
$$\mathbf E=\frac{q}{4\pi\epsilon_0 |\mathbf r|^3}\hat{\mathbf r}$$
This can be derived easily from integral form of Gauss’s law. Taking $V$ as a sphere of radius $r$ centered at the point charge,
$$\int_{\partial V}\mathbf E\cdot d\mathbf A=\frac{q}{\epsilon_0}$$
$$4\pi r^2 E=\frac{q}{\epsilon_0}\implies E=\frac{q}{4\pi\epsilon_0 r^2}$$

However, I can’t see how this can be derived directly from the differential form of Maxwell’s equations. Maxwell’s equations only says the divergence and curl of electric field is zero everywhere away from the point charge, which isn’t quite useful at all.

How can the formula for electric field of a point charge derived directly from differential form of Maxwell’s equations? 
(Note: it is of course possible to convert the differential form into integral form and the proceed as I’ve shown. However, that’s somewhat trivial and I do not consider that as a direct derivation.)
 A: When you have a point charge, the charge density is not really defined, as it is basically infinite at the origin, and zero everywhere else. You would need a distribution rather than a function to describe it. The derivation here is not completely rigorous as I am treating $\delta_{\mathbb{R}^3}(\vec{r})$ as a normal function when doing the integrations...
We can use the Dirac delta function over $\mathbb{R}^3$ to write $$\rho(\vec{r}) = q\delta_{\mathbb{R}^3}(\vec{r})$$
Then the Gauss equation gives $$\nabla \cdot \vec{E} = \frac{q}{\varepsilon_0}\delta_{\mathbb{R}^3}(\vec{r})$$
Since the problem is spherically symmetric, we compute the divergence in spherical coordinates and use $\vec{E} = E_r \hat{r}$ (where $E_r$ is the radial component of $\vec{E}$ and $\hat{r}$ is a unit radial vector): $$\nabla \cdot \vec{E} = \frac{1}{r^2}\frac{d}{dr}(r^2E_r)$$
So we get Gauss's equation as:
$$\frac{d}{dr}(r^2E_r) = \frac{q}{\varepsilon_0}r^2\delta_{\mathbb{R}^3}(r\hat{r}) \tag{*}$$
For $r \gt 0$ the equation is just $\frac{d}{dr}(r^2E_r) = 0$, so that $r^2E_r = C$ for some constant $C$.
Therefore $E_r = \frac{C}{r^2}$ for $r \gt 0$.
To find the actual value of $C$, we integrate (*) from $0$ to $\infty$ and use $r^2E_r = C$ to find $$C = \frac{q}{\varepsilon_0}\int_{0}^{\infty}\delta_{\mathbb{R}^3}(r\hat{r})r^2dr$$ (note that this is not fully rigorous since the usual Fundamental theorem of calculus applies to functions, not distributions).
Now the main property of $\delta_{\mathbb{R}^3}(\vec{r})$ is that $$\iiint\limits_{\mathbb{R}^3} \delta_{\mathbb{R}^3}(\vec{r}) dV = 1$$
Using spherical coordinates, this is $$1 = \int\limits_{0}^{\infty}\int\limits_{-\pi}^{\pi}\int\limits_{0}^{2\pi}\delta_{\mathbb{R}^3}(r\hat{r})r^2\sin\varphi \text{ d}\theta\text{ d}\varphi \text{ d}r = \int\limits_{0}^{\infty}\delta_{\mathbb{R}^3}(r\hat{r})r^2\text{ d}r \int\limits_{-\pi}^{\pi}\sin\varphi\text{ d}\varphi\int\limits_{0}^{2\pi}d\theta = 4\pi\int\limits_{0}^{\infty}\delta_{\mathbb{R}^3}(r\hat{r})r^2\text{ d}r$$
Therefore $C = \frac{q}{4\pi\varepsilon_0}$ and we obtain $$\vec{E} = \frac{q}{4\pi\varepsilon_0r^2}\hat{r}$$
On the other hand, this derivation is not actually too different from the one using the integral version of Gauss's theorem, since we performed an integration and used the 1-d fundamental theorem of calculus to find $C$, which is equivalent to using the divergence theorem in 3-d.
A: It'd be hard to get the multiplicative factor, I think, but the basic principle is not hard:
The divergence in spherical coordinates is handily provided by Wikipedia as$$\nabla\cdot A = 
{1 \over r^2}{\partial \left( r^2 A_r \right) \over \partial r}
+ {1 \over r\sin\theta}{\partial \over \partial \theta} \left(  A_\theta\sin\theta \right)
+ {1 \over r\sin\theta}{\partial A_\varphi \over \partial \varphi}.
$$Looking for a spherically symmetric field gives just $$
{d \left( r^2 E_r \right) \over d r} = r^2~ {\rho(r)\over\epsilon_0}.
$$
To get the proportionality right one would want to solve for $\rho(r)=\{3q/(4\pi R^3)\text{ if } r< R\text{ else }0\}$ with boundary condition $E_r(r=0)=0.$
