Solution (source code)

= Solution

Let $h=\Delta x$. A one-dimensional <Taylor expansion> gives
$$
u(x+h)-2u(x)+u(x-h)
=h^2u_{xx}(x)+\frac{h^4}{12}u_{xxxx}(x)+O(h^6).
$$
Adding the three coordinate directions shows that the exact solution has stencil defect
$$
L_hu-h^2f
=\frac{h^4}{12}(u_{x_1x_1x_1x_1}+u_{x_2x_2x_2x_2}+u_{x_3x_3x_3x_3})+O(h^6)
=O(h^4).
$$
This is the local defect of the <seven-point Dirichlet Laplacian>. Its inverse has max-norm size $O(h^{-2})$: this follows from the <discrete maximum principle>, for example by comparison with the grid restriction of $x_1(1-x_1)+x_2(1-x_2)+x_3(1-x_3)$, whose stencil is $-6h^2$. Therefore the grid error is $O(h^2)$.

The question writes this error as $O(h^{p+1})$, so
$$
\boxed{p+1=2,\qquad p=1.}
$$
In the more usual terminology, the <finite difference method> is second-order accurate.