  Allowing for space- and time-dependence of mass in Klein--Gordon equations
resolves the problem of negative probability density and of violation of Lorenz
covariance of interaction in quantum mechanics. Moreover it extends their
applicability to the domain of quantum cosmology, where the variation in mass
may be accompanied by high oscillations. In this paper we propose a third-order
exponential integrator, where the main idea lies in embedding the oscillations
triggered by the possibly highly oscillatory component intrinsically into the
numerical discretisation. While typically high oscillation requires
appropriately small time steps, an application of Filon methods allows
implementation with large time steps even in the presence of very high
oscillation. This greatly improves the efficiency of the time-stepping
algorithm.
  Proof of the convergence and its rate are nontrivial and require alternative
representation of the equation under consideration. We derive careful bounds on
the growth of global error in time discretisation and prove that, contrary to
standard intuition, the error of time integration does not grow once the
frequency of oscillations increases. Several numerical simulations are
presented to confirm the theoretical investigations and the robustness of the
method in all oscillatory regimes.
