The Galerkin method offers a powerful tool in the solution of differential equations and function approximation on the real interval [Γ1, 1]. By expanding the unknown function in appropriately chosen global basis functions, each of which explicitly satisfies the given boundary conditions, in general this scheme converges exponentially fast and almost always supplies the most terse representation of a smooth solution. To date, typical schemes have been defined in terms of a linear combination of two Jacobi polynomials. However, the resulting functions do not inherit the expedient properties of the Jacobi polynomials themselves and the basis set will not only be non-orthogonal but may, in fact, be poorly conditioned. Using a Gram-Schmidt procedure, it is possible to construct, in an incremental fashion, polynomial basis sets that not only satisfy any linear homogeneous boundary conditions but are also orthogonal with respect to the general weighting function Γ°1 Γ xΓ a Γ°1 ΓΎ xΓ b . However, as it stands, this method is not only cumbersome but does not provide the structure for general index n of the functions and obscures their dependence on the parameters Γ°a; bΓ. In this paper, it is shown that each of these Galerkin basis functions, as calculated by the Gram-Schmidt procedure, may be written as a linear combination of a small number of Jacobi polynomials with coefficients that can be determined. Moreover, this terse analytic representation reveals that, for large index, the basis functions behave asymptotically like the single Jacobi polynomial P Γ°a;bΓ n Γ°xΓ. This new result shows that such Galerkin bases not only retain exponential convergence but expedient function-fitting properties too, in much the same way as the Jacobi polynomials themselves. This powerful methodology of constructing Galerkin basis sets is illustrated by many examples, and it is shown how the results extend to polar geometries. In exploring more generalised definitions of orthogonality involving derivatives, we discuss how a large class of differential operators may be discretised by Galerkin schemes and represented in a sparse fashion by the inverse of band-limited matrices.