In this paper, we study an accurate and stable numerical method for the system of partial differential equations defined by the linearized Westervelt equation and the Pennes bioheat equation. This system was proposed in the literature to describe the high-intensity focused ultrasound (HIFU) propagation through a target tissue and the correspondent temperature evolution due to the conversion of acoustic energy into heat. In the bioheat equation, the heat source term is allowed to depend on the acoustic field through \( p^2, p^2_t \) or on suitable time-averaged versions of these quantities. The method can be seen simultaneously as a finite difference method and as a fully discrete in space piecewise linear finite element method. We establish second order convergence in space for the discrete time derivative and for the discrete gradient of the numerical solution and first order convergence in time, including the case where the temperature equation is driven by the aforementioned heat sources. Semi-discrete and fully discrete in time and space discretizations are studied using the Bramble-Hilbert lemma to avoid the usual smoothness requirements on the solution. Simulation results illustrating the convergence estimates and the qualitative behavior of the solution of the differential problem are also included.