Obtaining PDE backstepping controller or observer gains requires the solution of kernel PDEs - one or more hyperbolic PDEs on a triangular (Goursat) domain with non-standard boundary conditions. The numerical solution of these equations is a challenge that every designer applying backstepping must eventually face, except for the simplest of cases where explicit solutions are available. In addition, recent backstepping designs for coupled systems exhibit discontinuous behavior which must be accurately captured with the numerical approximation. In this paper, we propose a power series method as an alternative to other approaches. This method, which was introduced years ago in combination with convex optimization but whose convergence has only been recently formally established, offers several advantages, most of all simplicity, as it is quite easy to grasp and implement. Other features of the method include precision, adaptability to settings with discontinuous kernels, and the ability to produce symbolical kernels depending on parameters. The paper provides the necessary theoretical background, which borrows fundamental results from complex analysis and leverages already written kernel well-posedness proofs to show the existence of a power series solution. Links to the codes used in the simple examples are given; these are easily adaptable Mathematica notebooks. A complex multi-kernel, multiple-discontinuity example used in the stabilizing feedback law of a multilayer Timoshenko beam is given at the end, demonstrating the applicability of the method to some challenging families of kernel equations appearing in recent backstepping designs.