For the precision calculations in perturbative Quantum Chromodynamics (QCD) gigantic expressions (several GB in size) in terms of highly complicated divergent multi-loop Feynman integrals have to be calculated analytically to compact expressions in terms of special functions and constants. In this article we derive new symbolic tools to gain large-scale computer understanding in QCD. Here we exploit the fact that hypergeometric structures in single and multiscale Feynman integrals emerge in a wide class of topologies. Using integration-by-parts relations, associated master or scalar integrals have to be calculated. For this purpose it appears useful to devise an automated method which recognizes the respective (partial) differential equations related to the corresponding higher transcendental functions. We solve these equations through associated recursions of the expansion coefficient of the multivalued formal Taylor series. The expansion coefficients can be determined using either the package in the case of linear difference equations or by applying heuristic methods in the case of partial linear difference equations. In the present context a new type of sums occurs, the Hurwitz harmonic sums, and generalized versions of them. The code transforming classes of differential equations into analytic series expansions is described. Also partial difference equations having rational solutions and rational function solutions of Pochhammer symbols are considered, for which the code is designed. Generalized hypergeometric functions, Appell-, Kampé de Fériet-, Horn-, Lauricella-Saran-, Srivasta-, and Exton–type functions are considered. We illustrate the algorithms by examples.