The paper introduces k-Wave, a MATLAB toolbox designed for the simulation and reconstruction of photoacoustic wave fields. k-Wave is designed to simplify and expedite realistic photoacoustic modeling, particularly in heterogeneous media. The toolbox includes functions for forward simulations based on a k-space pseudo-spectral time domain solution to coupled first-order acoustic equations in one, two, and three dimensions. It also supports time reversal image reconstruction for arbitrarily shaped measurement surfaces and a one-step image reconstruction algorithm for planar detector geometries using the fast Fourier transform (FFT). The paper discusses the architecture and usage of k-Wave, highlighting its capabilities in improving time reversal reconstructions with sparse detector arrays, comparing time reversal and FFT-based reconstruction methods, and optimizing computational speed through parallel execution using a graphics processing unit (GPU). The toolbox is freely available and can be used for various applications in acoustics and ultrasonics, including conventional diagnostic ultrasound, seismology, and environmental noise propagation.The paper introduces k-Wave, a MATLAB toolbox designed for the simulation and reconstruction of photoacoustic wave fields. k-Wave is designed to simplify and expedite realistic photoacoustic modeling, particularly in heterogeneous media. The toolbox includes functions for forward simulations based on a k-space pseudo-spectral time domain solution to coupled first-order acoustic equations in one, two, and three dimensions. It also supports time reversal image reconstruction for arbitrarily shaped measurement surfaces and a one-step image reconstruction algorithm for planar detector geometries using the fast Fourier transform (FFT). The paper discusses the architecture and usage of k-Wave, highlighting its capabilities in improving time reversal reconstructions with sparse detector arrays, comparing time reversal and FFT-based reconstruction methods, and optimizing computational speed through parallel execution using a graphics processing unit (GPU). The toolbox is freely available and can be used for various applications in acoustics and ultrasonics, including conventional diagnostic ultrasound, seismology, and environmental noise propagation.