Rendering highly scattering participating media using brute force path tracing is a challenge. The diffusion approximation reduces the problem to solving a simple linear partial differential equation. Flux-limited diffusion introduces nonlinearities to alleviate the approximation error but introduces several ad-hoc assumptions. Both methods are based on the spherical harmonics expansion of the radiance field, that is truncated after the first order. In this paper, we investigate the open question of whether higher orders provide a viable alternative to these two approaches. Increasing the order introduces a set of increasingly complex coupled partial differential equations, whose growing number and complexity make them very difficult to work with. We use a computer algebra framework for representing and manipulating the underlying mathematical equations and use it to derive the time-independent real-valued P_N-equations for arbitrary orders. We further present a staggered-grid P_N-solver and generate its stencil code directly from the expression tree of the P_N-equations. Finally, we discuss how our method compares to prior work for various standard problems.