A method for the numerical evaluation of the swallowtail canonical integral S(x, y, z) and its partial derivatives ∂S/∂x, ∂S/∂y and ∂S/∂z is described. These integrals are required for applications in semiclassical collision theory and other short wavelength phenomena. The method integrates a set of three linear ordinary differential equations satisfied by S(x, y, z). We first explain the method for the simpler cases of the Airy and Pearcey integrals. For S(x, y, z), the integration starts at (0, 0, 0) and then proceeds in the x direction to (x, 0, 0). This is followed by an integration in the y direction to (x, y, 0). Finally integration in the z direction takes us to (x, y, z). The initial conditions at (0, 0, 0) are obtained from the exact series representation of S(x, y, z). The method has been implemented on a computer and a discussion of its advantages and disadvantages is given. Isometric plots of |S(x, y, z)| are presented for fixed values of x, y and z in order to illustrate in a systematic way the main properties of the swallowtail integral. A knowledge of the number of real stationary phase points in different regions of (x, y, z) space together with the shape of the swallowtail caustic are a useful way of understanding and rationalizing the structure in the |S(x, y, z)| plots.