This paper proposes an algorithm based on the Shultz iterative method and divided differences to find the roots of nonlinear equations. Using the new method, we present a rapid and powerful algorithm to compute an approximate inverse of an invertible tensor. Analysis of the convergence error shows that the convergence order of the method is a linear combination of the Fibonacci sequence and also is rapid and powerful in finding and keeping sparsity of the obtained approximate inverse of the sparse tensors. The algorithm is extended for computing the Moore-Penrose inverse of a tensor. As an application, we use the iterates obtained by the algorithm as a preconditioner for the tensorized Krylov subspace method, e.g., LSQR based tensor form to solve the multilinear system A N X = B. Several examples are also provided to show the efficiency of the proposed method. Finally, some concluding remarks are given.