I think the code should be organised as follows:

sympy/functions/special/bessel.py
To this file I will add classes representing the Bessel type function families. (The file aready implements spherical integral order bessel functions, which are expressed directly using trigonometric functions.)

sympy/functions/special/hyper.py
This new file will contain classes for the Hypergeometric and Meijer G function families. The method _eval_expand_func() will call hyperexpand() (see below) to (try to) write a Hypergeometric or Meijer G function in terms of named special functions. Hypergeometric and Meijer G functions will not expand themselves automatically. The class for Meijer G functions probably also deserves a few helper methods to aid implementation of the integration algorithm (e.g. a method to absorb factors). I will implement these as need arises. The Hypergeometric function class will also implement some special values (i.e. , etc for sets of parameters where closed formulae exist). Only the simplest cases will be treated here, since there is vast literature on such “Hypergeometric identities”. If need ever arises, Zeilberger’s algorithm can be used to find many new identities.

sympy/simplify/hyperexpand.py
This new file will contain right at the top a “knowledge base” of Hypergeometric (and eventually Meijer G) function identities. This should be in particularly readable form, so that new identities can be added easily. The rest of the file will process this knowledgebase to create the data structures necessary for the function expansion algorithm (at module load time), and actually implement the hyperexpand() function that does the rewriting of Hypergeometric and Meijer G functions in terms of known special functions.

sympy/integrals/meijer.py
This new file will implement a function meijerint() that tries to do definite integration using the Meijer G function transform. An important part of this is the rewriting of essentially arbitrary functions in terms of Meijer G functions. I think this should be in this file, since it essentially amounts to computing Mellin transforms, i.e. is tightly coupled to integration. The rewrite code will reuse the knowledge base from hyperexpand.py.

sympy/integrals/integrals.py
The function integrate() will be changed to call meijerint() for definite integrals, if everything else failed.
I will also implement randomised unit tests, and add numerical tests for all function transformations and knowledge base entries.
Later on more classes of special functions will be added, and the knowledge base will be extended.
I intend to create my first pull request after the addition of the Bessel, Hypergeometric and Meijer G function classes (week one), and the second one after the expansion of Hypergeometric functions works reasonably well (hopefully after week three – but note that there will be my oneweek break in this time so the week numbers will shift). The third pull request should include working integration code (hopefully after week four). Planning ahead more seems hard.
The functions in bessel.py could use some cleaning up. First off, they should be implemented nonrecursively (ideally using the Polys, like we do for Legendre polynomials). Second off, we need actual classes, so we can represent the coefficients with symbolic indices (we could also do fun stuff with this like add solvers to ode.py that return the respective coefficient function given the defining ODE).
“Hypergeometric and Meijer G functions will not expand themselves automatically.” Good. This is almost always the best strategy.
We’ll see how the performance of meijerint() is, but likely it should be called before heurisch() (but after the other ones), because heurisch() is very slow.
For randomized tests, see issue 2281.
And don’t worry about planning ahead too much 🙂