The numerical approximation of the Laplace equation with inhomogeneous mixed boundary conditions in 2D with lowest-order Raviart-Thomas mixed finite elements is realized in three flexible and short MATLAB programs. It is the aim of this paper to derive, document, illustrate, and validate the three MATLAB implementations EBmfem, LMmfem, and CRmfem for further use and modification in education and research. A posteriori error control with a reliable and efficient averaging technique is included to monitor the discretization error. Therein, emphasis is on the correct treatment of mixed boundary conditions. Numerical examples illustrate some applications of the provided software and the quality of the error estimation.