Following Zeilberger's holonomic systems approach, it was F. Chyzak who laid the algorithmic foundations for the automatic manipulation of multivariate "holonomic" (\partial-finite) functions in his PhD thesis in 1998. We believe that his work did not yet receive the attention that it deserves. In the past year we have implemented his and related algorithms in Mathematica in the package HolonomicFunctions. In this talk, we shall present several extensions of Chyzak's work that turned out to be very useful. The main focus will be on recent results which were obtained in cooperation with J. Schoeberl (RWTH Aachen) and which play a key role in the numeric simulation of electromagnetic waves (using high order discontinuous Galerkin finite lements). Our contribution was to find certain transformation formulae for the basis functions which are basically products of Legendre and Jacobi polynomials.