The p-Laplacian system for the electromagnetic effects in a high-temperature superconductor (HTS) is known to be tricky to solve because of its nonlinearity and stiffness which are caused by the power-law relation in the density function. The stiffness leads the parabolic system to have solutions which are almost discontinuous like shocks in hyperbolic problems. Noticing this fact, Giovanni Naldi modified the relaxation method of Shi Jin for porous media problem and showed improvements in accuracy and larger time-steps. However, applying the method of Naldi requires the p-Laplacian system to be differentiated first and then recovered by integration, and this integration has massive truncation errors when p is large. Therefore it is rather better to use the original relaxation method of Shi Jin to solve p-Laplcian problem. However the method needs very small mesh size in time and shows some diffusions. To resolve this problem we modified the relaxation system using the underlying idea of Naldi's method. The new method is expected to give slower wave speeds which would enlarge the mesh size in time, and have a much more suitable form than Naldi’s method for stiff problems when p is large.