This paper aims to develop efficient numerical methods for computing the inverse of matrix \(\varvec{\varphi }\) -functions, \(\varvec{\psi }_{\varvec{\ell }}\varvec{(A)}\varvec{:=}\varvec{ (\varphi }_{\varvec{\ell }}\varvec{(A))}^{\varvec{-1}}\) , for \(\varvec{\ell =1,2,\ldots ,}\) when \(\varvec{A}\) is a large and sparse matrix with eigenvalues in the open left half-plane. While \(\varvec{\varphi }\) -functions play a crucial role in the analysis and implementation of exponential integrators, their inverses arise in solving certain direct and inverse differential problems with non-local boundary conditions. We propose an adaptation of the standard scaling-and-squaring technique for computing \(\varvec{\psi }_{\varvec{\ell }}\varvec{(A)}\) , based on the Newton-Schulz iteration for matrix inversion. The convergence of this method is analyzed both theoretically and numerically. In addition, we derive and analyze Padé approximants for approximating \(\varvec{\psi }_{\varvec{1}}\varvec{(A/2}^{\varvec{s}}\varvec{)}\) , where s is a suitably chosen integer, necessary at the root of the squaring process. Numerical experiments demonstrate the effectiveness of the proposed approach.
