Laplace operator:
$\Delta f = \nabla \cdot \nabla f = 4 \pi \rho$
Let $f = \frac{c}{|x|}$
So $\nabla_a f = -c \frac{1}{|x|^2} \frac{x_a}{|x|}$
So $\nabla^a \nabla_a f = -c g^{ab} ( \nabla_b (\frac{1}{|x|^2}) \frac{x_a}{|x|} + \frac{1}{|x|^2} \nabla_b \frac{x_a}{|x|} )$
$= -c g^{ab} ( -2 \frac{1}{|x|^3} \frac{x_b}{|x|} \frac{x_a}{|x|} + \frac{1}{|x|^2} ( g_{ab} \frac{1}{|x|} - x_a \frac{1}{|x|^2} \frac{x_b}{|x|} ) )$
$= -c g^{ab} ( \frac{1}{|x|^3} g_{ab} - \frac{3}{|x|^5} x_a x_b )$
$= -c (n-3) \frac{1}{|x|^3}$ for $n =$ the dimension of the manifold

scalar fields in Cartesian coordinates:

$\Delta f = \nabla \cdot \nabla f = \Sigma_i \frac{\partial^2 f}{\partial (x^i)^2}$

finite-difference:
$\Delta f \approx \Sigma_i \frac{ (\frac{f(x + e_i) - f(x)}{\Delta x_i}) - (\frac{f(x) - f(x - e_i)}{\Delta x_i}) }{\Delta x_i}$
$\approx \Sigma_i \frac{ f(x + e_i) - 2 f(x) + f(x - e_i) }{(\Delta x_i)^2}$
$\approx \Sigma_i \frac{ f(x + e_i) + f(x - e_i) }{(\Delta x_i)^2} - 2 f(x) \Sigma_i \frac{1}{(\Delta x_i)^2}$
$\approx \Sigma_i \frac{ f_{I(x + e_i)} + f_{I(x - e_i)} }{(\Delta x_i)^2} - 2 f_{I(x)} \Sigma_i \frac{1}{(\Delta x_i)^2}$
$\approx A_{jk} f_j$
for $A_{jj} = -2 \Sigma_i \frac{1}{(\Delta x_i)^2}$
for $A_{jk} = \frac{1}{(\Delta x_i)^2}$ for $j = I(x), k = I(x + e_i)$, so notice that no summation over $i$ is used.

scalar fields in curvilinear coordinates:

$\Delta f = \nabla \cdot \nabla f = f_{;ab} g^{ab}$
$= (f_{,ab} - {\Gamma^c}_{ab} f_{,c}) g^{ab}$
$= f_{,ab} g^{ab} - \Gamma^c f_{,c}$

another derivation:
$\Delta \phi = g^{ab} \nabla_a \nabla_b \phi$
$= g^{ab} \nabla_b e_a(\phi)$
$= g^{ab} (e_b( e_a(\phi)) - {\Gamma^c}_{ba} e_c(\phi))$
$= g^{ab} e_b(e_a(\phi)) - {\Gamma^{ab}}_b e_a(\phi)$
$= {\Gamma_b}^{ab} e_a(\phi) - {\Gamma_b}^{ab} e_a(\phi) - {\Gamma^{ab}}_b e_a(\phi) + g^{ab} e_a(e_b(\phi)) $
$= {\Gamma_{cb}}^c g^{ab} e_a(\phi) - g^{bd} \Gamma_{dcb} g^{ac} e_a(\phi) - g^{bd} \Gamma_{cbd} g^{ac} e_a(\phi) + g^{ab} e_a(e_b(\phi)) $
$= {\Gamma^c}_{bc} g^{ab} e_a(\phi) - g^{ac} (\Gamma_{dcb} + \Gamma_{cbd}) g^{db} e_a(\phi) + g^{ab} e_a(e_b(\phi)) $
$= {\Gamma^c}_{bc} g^{ab} e_a(\phi) - g^{ac} ( \Gamma_{dbc} + c_{cbd} + \Gamma_{cbd} ) g^{db} e_a(\phi) + g^{ab} e_a(e_b(\phi)) $
$= \frac{1}{\sqrt{|g|}} e_b(\sqrt{|g|}) g^{ab} e_a(\phi) - g^{ac} (e_b(g_{cd}) + c_{cbd}) g^{db} e_a(\phi) + g^{ab} e_a(e_b(\phi)) $
$= \frac{1}{\sqrt{|g|}} ( e_b(\sqrt{|g|}) g^{ab} e_a(\phi) + \sqrt{|g|} e_b(g^{ab}) e_a(\phi) + \sqrt{|g|} g^{ab} e_a(e_b(\phi)) - {c^{ab}}_b e_a(\phi) )$
$= \frac{1}{\sqrt{|g|}} e_b(\sqrt{|g|} g^{ab} e_a(\phi)) - \frac{1}{\sqrt{|g|}} {c^{ab}}_b e_a(\phi) $
TODO double-check the commutation code, or see if it can disappear somewhere.