  Elliptic partial differential equations (ePDEs) appear in a wide variety of
areas of mathematics, physics and engineering. Typically, ePDEs must be solved
numerically, which sets an ever growing demand for efficient and highly
parallel algorithms to tackle their computational solution. The Scheduled
Relaxation Jacobi (SRJ) is a promising class of methods, atypical for combining
simplicity and efficiency, that has been recently introduced for solving linear
Poisson-like ePDEs. The SRJ methodology relies on computing the appropriate
parameters of a multilevel approach with the goal of minimizing the number of
iterations needed to cut down the residuals below specified tolerances. The
efficiency in the reduction of the residual increases with the number of levels
employed in the algorithm. Applying the original methodology to compute the
algorithm parameters with more than 5 levels notably hinders obtaining optimal
SRJ schemes, as the mixed (non-linear) algebraic-differential equations from
which they result become notably stiff. Here we present a new methodology for
obtaining the parameters of SRJ schemes that overcomes the limitations of the
original algorithm and provide parameters for SRJ schemes with up to 15 levels
and resolutions of up to $2^{15}$ points per dimension, allowing for
acceleration factors larger than several hundreds with respect to the Jacobi
method for typical resolutions and, in some high resolution cases, close to
1000. Furthermore, we extend the original algorithm to apply it to certain
systems of non-linear ePDEs.
