A new numerical method for tackling the three-dimensional Heston–Hull–White partial differential equation (PDE) is proposed. This PDE has an application in pricing options when not only the asset price and the volatility but also the risk-free rate of interest are coming from stochastic nature. To solve this time-dependent three-dimensional PDE as efficiently as possible, high order adaptive finite difference (FD) methods are applied for the application of method of lines. It is derived that the new estimates have fourth order of convergence on non-uniform grids. In addition, it is proved that the overall procedure is conditionally time-stable. The results are upheld via several numerical tests.