Acoustic and elastic wave propagation are important phenomena in various engineering applications. An example in biomedical engineering that interests us is modeling the propagation of ultrasound fields in the human body. Treatments using focused ultrasound deliver acoustic energy to a specific region, causing thermal or mechanical effects. Computational simulations are essential in understanding the impact of bone on focal aberrations and local overheating. Simulating therapeutic ultrasound requires solving the Helmholtz equation for acoustics and the Navier equation for elastodynamics in unbounded regions at high frequencies. The dominant wavelengths at MHz frequencies are much smaller than the relevant regions of the skull and the rib cage that affect the ultrasound beam. Hence, one needs fine meshes and, therefore, fast numerical solvers. Boundary Element Methods (BEM) naturally handle unbounded domains and tend to scale well with higher frequencies. However, they must be coupled with volume discretizations to model locally heterogeneous domains. Coupling with the Finite Element Method (FEM) requires dedicated preconditioners to solve the linear system and requires fine meshes at high frequencies due to the pollution effect. We designed the Volume-Surface Integral Equations (VSIE) algorithm, which discretizes a set of coupled volume and surface integral equations. The volumetric equations handle heterogeneous material parameters. Their coupling with surface integrals, similar to the ones in BEM, handles unbounded domains and high-contrast interfaces. We validated the VSIE algorithm against analytical solutions, the BEM, and FEM-BEM approaches. The VSIE shows high numerical accuracy at relatively coarse meshes and stability at resonances. The algorithm is also robust with increasing frequencies and high-contrast interfaces between materials. Furthermore, hierarchical matrix compression stores the dense system matrix using a fraction of the memory required for the whole matrix. Finally, we will demonstrate that the VSIE can model realistic therapeutic ultrasound scenarios.