The p-Laplacian is a non-linear generalization of the Laplace operator. A discrete version of p-Laplacian, the graph p-Laplacian, has been successfully used in various applications, including data clustering, dimensionality reduction and other tasks, since non-linearity better captures the underlying geometry of the data. In this work we present a novel iterative computational approach to find the p-Laplacian eigenpairs which extends the inverse shifted power method, a standard numerical technique to find specific eigenvalues/eigenfunctions to the nonlinear p-Laplacian operator. This involves solving a related nonlinear problem at each step using quasi-Newton method with Broyden’s update. By varying the parameter σ in a suitably defined range, the method converges to the spectrum of the p-Laplacian, for 1 < p < +∞. This reveals a continuous deformation of the classical 2-Laplacian eigenfunctions, with modified sharpness, but also gives rise to novel eigenfunctions with distinct geometric profiles. The accuracy of the results is measured via Orthogonal Least Square method, which also returns the estimated eigenvalues.