  Generating functions of BPS invariants for N=4 U(r) gauge theory on a
Hirzebruch surface with r=2 and 3 are computed. The BPS invariants provide the
Betti numbers of moduli spaces of semi-stable 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.
