Simulating the propagation of elastic waves in multi-layered media has many applications. A common approach is to use matrix methods where the elastic wave-field within each material layer is represented by a sum of partial-waves along with boundary conditions imposed at each interface. While these methods are well-known, coding the required matrix formation, inversion, and analysis for general multi-layered systems is non-trivial and time-consuming. Here, a new open-source toolbox called ElasticMatrix is described which solves the problem of acoustic and elastic wave propagation in multi-layered media for isotropic and transverse-isotropic materials where the wave propagation occurs in a material plane of symmetry. The toolbox is implemented in MATLAB using an object oriented programming framework and is designed to be easy to use and extend. Methods are provided for calculating and plotting dispersion curves, displacement and stress fields, reflection and transmission coefficients, and slowness profiles.Keywords partial-wave method · global matrix method · elastic waves 1 Motivation and Significance Matrix models of wave propagation in multi-layered elastic solids have had a significant contribution to research areas such as acoustics, geophysics and electromagnetics. A few examples include: structural health monitoring [1], characterisation of interface bonding [2], detection of debonding in joints [3], measuring material properties [4], designing composite layered structures [5], mode sorting of guided waves [6], the physical interpretation of guided wave structures [7], modelling the directional response of Fabry-Pérot ultrasound sensors [8], reflection and transmission of plane waves [9], elastography of layered soft tissues [10], and ice detection on wind turbines [11].Matrix methods, in particular the partial-wave and global matrix method, represent the stress and displacement fields as a sum of partial-waves for each material of the layered-structure. Each partial-wave represents an upward or downward travelling (quasi-)compressional or (quasi-)shear wave. By invoking boundary conditions at the interfaces of adjacent layers, the partial-wave amplitudes and field properties of the first layer can be related to the last in the form of a 'global' matrix. The resulting matrix equation can be used in two different ways. Firstly, the roots of the equation can be found which give the modal solutions or dispersion curves. Secondly, a subset of partial-wave amplitudes can be defined and the remaining amplitudes solved for. This can be used to calculate the displacement and stress fields within the multi-layered structure when a plane wave is incident. This method will be discussed further in Section 2.Despite its usefulness, there are few available implementations of the partial-wave method. The current state-of-theart implementation is Disperse [12]. This software has been in development since 1990 and is primarily focused on calculating the dispersion solutions for multi-layered structures. The Disperse softwar...