A standard artificial compression (AC) method for incompressible flow is u(n+1)(epsilon) - u(n)(epsilon)/k + u(n+1)(epsilon) . del u(n+1)(epsilon) + 1/2u(n+1)(epsilon) del . u(n+1)(epsilon) +del p(n+1)(epsilon) - nu Delta u(n+1)(epsilon) = f, epsilon p(n+1)(epsilon) - p(n)(epsilon)/k + del . u(n+1)(epsilon) = 0 for, typically, epsilon=k (timestep). It is fast, efficient and stable with accuracy O(epsilon+k). For adaptive (and thus variable) timestep k(n) (and thus epsilon = epsilon(n)) its long time stability is unknown. For variable k, epsilon this report shows how to adapt a standard AC method to recover a provably stable method. For the associated continuum AC model, we prove convergence of the epsilon = epsilon(t) artificial compression model to a weak solution of the incompressible Navier-Stokes equations as epsilon = epsilon(t) -> 0. The analysis is based on space-time Strichartz estimates for a non-autonomous acoustic equation. Variable epsilon, k numerical tests in 2d and 3d are given for the new AC method.