Ren­der­ing highly scat­ter­ing par­ti­cip­at­ing me­dia us­ing brute force path tra­cing is a chal­lenge. The dif­fu­sion ap­prox­im­a­tion re­duces the prob­lem to solv­ing a simple lin­ear par­tial dif­fer­en­tial equa­tion. Flux-lim­ited dif­fu­sion in­tro­duces non­lin­ear­it­ies to al­le­vi­ate the ap­prox­im­a­tion er­ror but in­tro­duces sev­er­al ad-hoc as­sump­tions. Both meth­ods are based on the spher­ic­al har­mon­ics ex­pan­sion of the ra­di­ance field, that is trun­cated after the first or­der. In this pa­per, we in­vest­ig­ate the open ques­tion of wheth­er high­er or­ders provide a vi­able al­tern­at­ive to these two ap­proaches. In­creas­ing the or­der in­tro­duces a set of in­creas­ingly com­plex coupled par­tial dif­fer­en­tial equa­tions, whose grow­ing num­ber and com­plex­ity make them very dif­fi­cult to work with. We use a com­puter al­gebra frame­work for rep­res­ent­ing and ma­nip­u­lat­ing the un­der­ly­ing math­em­at­ic­al equa­tions and use it to de­rive the time-in­de­pend­ent real-val­ued P_N-equa­tions for ar­bit­rary or­ders. We fur­ther present a staggered-grid P_N-solv­er and gen­er­ate its sten­cil code dir­ectly from the ex­pres­sion tree of the P_N-equa­tions. Fi­nally, we dis­cuss how our meth­od com­pares to pri­or work for vari­ous stand­ard prob­lems.