If $a,b$ are positive or negative at the same time, then scaling will do the trick. Let $$x_1' = c x_1, x_2' = d x_2,$$ then
$$
\frac{\partial^2 u}{\partial x_1'^2} = \frac{1}{c^2}\frac{\partial^2 u}{\partial x_1^2}, \quad \frac{\partial^2 u}{\partial x_2'^2} = \frac{1}{d^2}\frac{\partial^2 u}{\partial x_2^2}.
$$
Comparing with the equation we have $c = 1/\sqrt{a}$ and $d = 1/\sqrt{b}$, the fundamental solution to
$$
a\frac{\partial^2 u}{\partial x_1^2}+ b\frac{\partial^2 u}{\partial x_2^2}= \frac{\partial^2 u}{\partial x_1'^2} + \frac{\partial^2 u}{\partial x_2'^2}= \delta_{(0,0)}(x_1,x_2)
$$
is
$$
E = \frac{1}{4\pi} \ln(x_1'^2 + x_2'^2) = \frac{1}{4\pi} \ln\left(\frac{x_1^2}{a} + \frac{x_2^2}{b}\right).
$$
To verify it, notice your operator $\mathcal{D}$ acting $u$ is
$$
\mathcal{D}u = \nabla \cdot (A\nabla u) = a\frac{\partial^2 u}{\partial x_1^2}+ b\frac{\partial^2 u}{\partial x_2^2}, \quad \text{ where } A = \begin{pmatrix}a &0\\0&b\end{pmatrix}.
$$
$\nabla$ is just $D$, and $\nabla \cdot $ is the divergence operator.
You can find that
$$
\nabla u = \frac{1}{4\pi}\left(\frac{2 b x_1}{b x_1^2+a x_2^2}, \frac{2 ax_2}{b x_1^2+a x_2^2}\right),
$$
and acting $A$ on $\nabla u$ we have:
$$
A\nabla u = \frac{1}{4\pi}\left(\frac{2 ab x_1}{b x_1^2+a x_2^2}, \frac{2 a bx_2}{b x_1^2+a x_2^2}\right).
$$
When $x_1,x_2\neq 0$, take divergence:
$$
\nabla \cdot (A\nabla u) = \frac{\partial }{\partial x_1} \left(\frac{2 ab x_1}{b x_1^2+a x_2^2}\right)+ \frac{\partial }{\partial x_2} \left(\frac{2 a bx_2}{b x_1^2+a x_2^2}\right) = 0.
$$
When $x_1,x_2 = 0$, you can refer to my answer in this question: Green's theorem and flux
and the divergence is a Dirac delta (for we have a $1/(4\pi)$ factor).
If $a,b$ have different signs, then this is not an elliptic operator anymore but 1D a wave operator, the solutions are plane wave will be rather given in D'Alembert formula, and the fundamental solution is given by the Heaviside step function $H(x_2 -\sqrt{-b/a} x_1 )$.