 Solution

Let $h=\Delta x$. A one-dimensional Taylor expansion givesAdding the three coordinate directions shows that the exact solution has stencil defectThis 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})$, soIn the more usual terminology, the finite difference method is second-order accurate.