The efficient inversion of matrix polynomials is a critical challenge in computational mathematics. We design a procedure to determine the inverse of matrices polynomial of multidimensional Laplace matrices. The method is based on eigenvector and eigenvalue expansions. The method is consistent with previously known expressions of the inverse discretized Laplacian in one spatial dimension (Vermolen et al., 2022). The formalism is further extended to obtain closed form expressions for time-dependent problems.