Generating functions of BPS invariants for N = 4 U (r) gauge theory on a Hirzebruch surface with r ≤ 3 are computed. The BPS invariants provide the Betti numbers of moduli spaces of semistable sheaves. The generating functions for r = 2 are expressed in terms of higher level Appell functions for a certain polarization of the surface. The level corresponds to the self-intersection of the base curve of the Hirzebruch surface. The non-holomorphic functions are determined, which added to the holomorphic generating functions provide functions, which transform as a modular form.