A new, freely available MATLAB toolbox for simulating and reconstructing photoacoustic wave fields is introduced. The toolbox, named k-Wave, is designed to simplify and accelerate realistic photoacoustic modeling. It uses a k-space pseudo-spectral time domain solution to coupled first-order acoustic equations for homogeneous or heterogeneous media in 1D, 2D, and 3D. The simulation functions can also be used as a flexible time reversal image reconstruction algorithm for arbitrary measurement surfaces. A one-step image reconstruction algorithm based on the fast Fourier transform (FFT) is also included. The architecture and use of the toolbox are described, along with several novel modeling examples. These include the use of data interpolation to improve time reversal reconstructions for sparse detector arrays, the comparison of time reversal and one-step FFT-based reconstruction algorithms, and performance enhancements via parallel execution using a graphics processing unit (GPU). The governing photoacoustic equations and their application to simulation and reconstruction are discussed. The k-space pseudo-spectral solution method and the use of a perfectly matched layer are also described. The k-Wave toolbox is introduced, and many of its functions are illustrated. Several novel simulation and reconstruction examples are presented, including the use of data interpolation for sparse detector arrays, a comparison of time reversal and one-step FFT-based image reconstruction algorithms, and performance enhancements via data casting and parallelization using the GPU. A summary and discussion of future work are given. The k-Wave toolbox is designed to make photoacoustic modeling easy and fast. The functions include simulation of photoacoustic wave fields, reconstruction of photoacoustic images, creation of geometric shapes, and utility functions. The simulation functions compute the time evolution of acoustic wave fields in homogeneous or heterogeneous media. The functions can also be used for time reversal image reconstruction. Additional image reconstruction functions allow the initial photoacoustic pressure to be estimated from data recorded over a linear or planar measurement surface. The geometry creation functions allow both Cartesian- and grid-based geometries to be defined. The utility functions are used for tasks such as grid creation, matrix smoothing, matrix interpolation, and file loading. The k-Wave toolbox is used for simulating and reconstructing photoacoustic wave fields. The simulation functions are based on a k-space solution to coupled acoustic equations. The time reversal image reconstruction functions are used to reconstruct the initial photoacoustic pressure distribution from recorded data. The one-step image reconstruction functions are used for planar measurement surfaces. The k-Wave toolbox is used for various applications, including the simulation of photoacoustic wave fields, the reconstruction of photoacoustic images, and the creation of geometric shapes. The toolbox is also used for optimizing computational speed through parallel execution using a GPU. The k-Wave toolbox is a valuable tool for researchers in the field of photoacoustic imaging.A new, freely available MATLAB toolbox for simulating and reconstructing photoacoustic wave fields is introduced. The toolbox, named k-Wave, is designed to simplify and accelerate realistic photoacoustic modeling. It uses a k-space pseudo-spectral time domain solution to coupled first-order acoustic equations for homogeneous or heterogeneous media in 1D, 2D, and 3D. The simulation functions can also be used as a flexible time reversal image reconstruction algorithm for arbitrary measurement surfaces. A one-step image reconstruction algorithm based on the fast Fourier transform (FFT) is also included. The architecture and use of the toolbox are described, along with several novel modeling examples. These include the use of data interpolation to improve time reversal reconstructions for sparse detector arrays, the comparison of time reversal and one-step FFT-based reconstruction algorithms, and performance enhancements via parallel execution using a graphics processing unit (GPU). The governing photoacoustic equations and their application to simulation and reconstruction are discussed. The k-space pseudo-spectral solution method and the use of a perfectly matched layer are also described. The k-Wave toolbox is introduced, and many of its functions are illustrated. Several novel simulation and reconstruction examples are presented, including the use of data interpolation for sparse detector arrays, a comparison of time reversal and one-step FFT-based image reconstruction algorithms, and performance enhancements via data casting and parallelization using the GPU. A summary and discussion of future work are given. The k-Wave toolbox is designed to make photoacoustic modeling easy and fast. The functions include simulation of photoacoustic wave fields, reconstruction of photoacoustic images, creation of geometric shapes, and utility functions. The simulation functions compute the time evolution of acoustic wave fields in homogeneous or heterogeneous media. The functions can also be used for time reversal image reconstruction. Additional image reconstruction functions allow the initial photoacoustic pressure to be estimated from data recorded over a linear or planar measurement surface. The geometry creation functions allow both Cartesian- and grid-based geometries to be defined. The utility functions are used for tasks such as grid creation, matrix smoothing, matrix interpolation, and file loading. The k-Wave toolbox is used for simulating and reconstructing photoacoustic wave fields. The simulation functions are based on a k-space solution to coupled acoustic equations. The time reversal image reconstruction functions are used to reconstruct the initial photoacoustic pressure distribution from recorded data. The one-step image reconstruction functions are used for planar measurement surfaces. The k-Wave toolbox is used for various applications, including the simulation of photoacoustic wave fields, the reconstruction of photoacoustic images, and the creation of geometric shapes. The toolbox is also used for optimizing computational speed through parallel execution using a GPU. The k-Wave toolbox is a valuable tool for researchers in the field of photoacoustic imaging.