Exact solutions for shear stress in a polymeric liquid subjected to large-amplitude oscillatory shear flow (LAOS) contain many Bessel functions. For the simplest of these, for instance, the corotational Maxwell fluid, in the closed form for its exact solution, Bessel functions appear 42 times, each of which is inside a summation. Approximate analytical solutions for shear stress in LAOS often take the form of the first few terms of a power series in the shear rate amplitude, and without any Bessel functions at all. There is thus practical interest in extending the Goddard integral expansion (GIE), to an arbitrary number of terms. In continuum theory, these truncated series are arrived at laboriously using the GIE. However, each term in the GIE requires much more work than its predecessor. For the corotational Maxwell fluid, for instance, the GIE for the shear stress has yet to be taken beyond the sixth power of the shear rate amplitude. In this paper, we begin with the exact solution for shear stress responses in corotational Maxwell fluids, and then perform an expansion by symbolic computation to confirm up to the sixth power, and to then continue the GIE. In this paper for example, we continue the GIE to the 40th power of the shear rate amplitude. We use Ewoldt grids to show our main result to be highly accurate. We also show the radius of convergence of the GIE to be infinite.