In this work, we present an adaptive Levin-type method for highprecision computation of highly oscillatory integrals with integrands of the form f (x) exp (i ωg(x)). If g has no real zero in the integration interval and the integrand is sufficiently smooth, the method can attain arbitrarily high asymptotic orders without computation of derivatives. Although the proposed method requires about 2d-digit working precision for a d-digit accurate approximation to the integral, it is more efficient than the numerical steepest decent method for some integrals. Moreover it produces sufficiently accurate approximations even if the integrand is only several times differentiable. For the implementation of the method we develop a Mathematica program, which can also be used for the well-known Levin-type method where the required derivatives are computed by Mathematica's automatic differentiation package. The effectiveness of the method is discussed in the light of a set of test integrals, one of which is computed to 1000 digit accuracy.