I don't know what you mean specifically by a "systematic way to do it." But to begin, the Dirac Delta and its derivatives are Generalized Functions, also known as Distributions.
Distributions are linear Functionals that map test functions from $C_C^\infty$ onto the reals. For the Dirac Delta, the functional definition is given as
$$\langle \delta_a, \phi \rangle =\phi(a) \tag 1$$
where $\phi$ is a suitable test function.
DERIVATIVE OF A DISTRIBUTION:
Suppose that $T$ is a distribution. We define the $n$'th derivative $T^{(n)}$ of $T$ such that for any $\phi\in C_C^\infty$, we have
$$\langle T^{(n)},\phi \rangle =(-1)^n\langle T,\phi^{(n)}\rangle \tag2$$
It is important to understand the difference between a function and its distributional counterpart. And it is, therefore, important to understand that differentiation of a function is not the same as differentiation of a distribution.
We can apply "systematic" rules such as the product rule for differention to distributions. Suppose that $T_1$ and $T_2$ are distributions and that the product $T_1T_2$ is a well-defined distribution (Note that the product of distributions might not be well-defined. For example, the product of the Dirac Delta and the Heaviside function is not defined.). Then, we have
$$\begin{align}
\langle (T_1T_2)',\phi\rangle&=-\langle T_1T_2,\phi'\rangle\\\\
&=-\langle T_1, T_2\phi'\rangle\\\\
&=-\langle T_1, (T_2\phi)'-T_2'\phi\rangle\\\\
&=\langle T_1',T_2\phi\rangle+\langle T_1,T_2'\phi\rangle\\\\
&=\langle T_1'T_2,\phi\rangle+\langle T_1T_2',\phi\rangle\\\\
&=\langle T_1T_2'+T_1T_2',\phi\rangle
\end{align}$$
Hence, in distribution, we assert that $(T_1T_2)'=T_1T_2'+T_1T_2'$.
In the following examples, we will systematically apply $(2)$ to determine the distributional interpretations of the (i) derivative of the Heaviside function $H$ and (ii) second mixed partials of the scalar potential function $T=\frac1{|\vec x|}$.
EXAMPLE $(1)$:
Let $H(t)$ denote the Heaviside (or unit step) function. The derivative of $H(t)$ is $0$ for $t\ne 0$ and fails to exist at $0$.
As a distribution, the derivative, $H'$ of $H$ can be evaluated as follows. Let $\phi \in C_C^\infty$. Then, applying $(2)$ we have
$$\begin{align}
\langle H',\phi\rangle&=-\langle H, \phi'\rangle\\\\
&=-\int_{-\infty}^\infty H(t)\phi'(t)\,dt\\\\
&=-\int_0^\infty \phi'(t)\,dt\\\\
&=\phi(0)\tag3
\end{align}$$
whence we see from $(3)$ that in distribution $H'=\delta_0$
EXAMPLE $2$:
As an example, define the function $T$ as $T(|\vec x|)=\frac1{|\vec x|}$ for $|\vec x|\ne 0$. The second mixed partial derivatives for $|\vec x|\ne 0$ are given by
$$\frac{\partial^2}{\partial x_i \partial x_j}\left(\frac1{|\vec x|}\right)=\frac{3x_ix_j-|\vec x|^2\delta_{ij}}{|\vec x|^5}$$
where $\delta_{ij}$ is the Kronecker Delta. Irrespective of how $T(0)$ is defined (or not defined), the second mixed partials of $T(|\vec x|)$ fail to exist at $0$.
Now, let $T(|\vec x|)$ be a distribution. Using $(2)$, we see that the second partial derivatives of $T$ are given by
$$\begin{align}
\langle \partial_{ij}T, \phi\rangle &=\langle T, \partial_{ij}\phi\rangle\\\\
&=\int_{\mathbb{R}^3}\frac1{|\vec x|} \frac{\partial^2 \phi}{\partial x_i \partial x_j}\,d^3\vec x\\\\
&=-\int_{\mathbb{R}^3}\frac{\partial}{\partial x_i}\left( \frac1{|\vec x|}\right)\frac{\partial \phi}{ \partial x_j}\,d^3\vec x\\\\
&=-\lim_{\varepsilon\to 0^+}\int_{\mathbb{R}^3\setminus B_\varepsilon(0)}\frac{\partial}{\partial x_i}\left( \frac1{|\vec x|} \right)\frac{\partial \phi}{ \partial x_j}\,d^3\vec x\tag4
\end{align}$$
where $B_\varepsilon(0)$ is a sphere of radius $\varepsilon$ centered at $0$.
Then integrating by parts the integral on the right-hand side of $(4)$ reveals
$$\begin{align}
\langle \partial_{ij}T, \phi\rangle &=-\phi(0)\int_0^{2\pi}\int_0^\pi n_i n_j \,\sin(\theta)d\theta\,d\phi+\lim_{\varepsilon\to 0^+}\int_{\mathbb{R}^3\setminus B_\varepsilon(0)}\frac{\partial^2}{\partial x_i\partial x_j}\left(\frac1{|\vec x|}\right)\phi\,d^3\vec x\\\\
&=-\frac{4\pi}{3}\phi(0)\delta_{ij}+\lim_{\varepsilon\to 0^+}\int_{\mathbb{R}^3\setminus B_\varepsilon(0)}\frac{3x_ix_j-|\vec x|^2\delta_{ij}}{|\vec x|^5}\phi\,d^3\vec x\tag5
\end{align}$$
where $n_i$ is the $i$'th component of the outer unit normal to $B_\varepsilon(0)$.
Therefore, we see from $(5)$ and $(1)$ that in distribution we have
$$\frac{\partial^2}{\partial x_i \partial x_j}\left(\frac1r\right)=-\frac{4\pi}{3}\delta_{ij}\delta_0(\vec x)+\text{PV}\left(\frac{3x_ix_j-|x|^2\delta_{ij}}{|\vec x|^5}\right)$$
SPECIAL CASES:
For $i=j$, we see that the second derivative of $\frac1{|\vec x|}$ with respect to $x$ is
$$\frac{\partial^2}{\partial x^2}\left(\frac1{|\vec x|}\right)=-\frac{4\pi}{3}\delta_0(\vec x)+\text{PV}\left(\frac{2x^2-y^2-z^2}{|\vec x|^5}\right)$$
And the Laplacian is $-4\pi \delta_0(\vec x)$ as expected.
where $C_{ij}$ is given by
$$C_{ij}=\int_0^{2\pi}\int_0^\pi n_i n_j ,\sin(\theta)d\theta,d\phi$$and we see that there is a Dirac Delta component along with a Principal Value term. So, I don't see a systematic way of evaluating a distributional derivative from methods that don't account for the distributions.
– Mark Viola Jun 11 '21 at 15:06