We propose a finite difference method to approximate weak distributional solutions of elliptic equations in the double-divergence form. Under minimal regularity assumptions on the coefficients, we resort to a regularisation argument. It turns our problem into a linear equation in the non-divergence form, with smooth coefficients. Regularity estimates build upon classical methods (e.g., Lax Equivalence Theorem) to yield convergence of the numerical method. To validate our strategy, we present three numerical examples in the planar setting. As far as we know, this is the first finite difference method for these equations, and our approach extends naturally to broader classes of models with low-regularity data.