The Marchenko method can be applied on cross-sections of 3D data using a 2D algorithm, but only a full 3D implementation can properly retrieve all 3D effects present in the data. The 3D implementation of the iterative Marchenko method is in principle a straightforward extension of the 2D method, it only requires an additional surface integration dimension. We present a 3D implementation of the Marchenko method based on an earlier implementation of the 2D Marchenko method. In the discussed implementation, special attention is given to an efficient kernel implementation and limiting the amount of data to read in, by using data compression. The algorithm properly handles 3D events and this is illustrated in the intermediate results of individual iterations of the method. Due to model-time and data-size constraints, the aperture of the data is limited, which leads to finite aperture artefacts that cause a lower convergence rate. We demonstrate the 3D method by retrieving the Green's functions for two models and comparing these functions to reference solutions. The two models are a horizontally layered medium and the complex SEG/EAGE overthrust model. Using the Marchenko method, imaging is applied to both models to show that false images caused by internal multiples are attenuated. Our 3D implementation is fully opensource and is intended to be used in future studies.