In this chapter, we present fv-unsat, a multipoint finite-volume–based solver for unsaturated flow in deformable and nondeformable porous media. The latter is described using the mixed form of Richards’ equation, whereas the former by the equations of unsaturated poroelasticity. The module aims at flexibility, relying heavily on discrete operators and equations, exploiting the automatic differentiation framework provided by the MATLAB Reservoir Simulation Toolbox (MRST). Our examples cover two numerical convergence tests and two three-dimensional practical applications, including the water infiltration process in a nondeformable soil column and a realistic desiccation process of a deformable clay sample using atmospheric boundary conditions. The resulting convergence rates are in agreement with previously reported rates for single-phase models, and the practical applications capture the physical processes accurately.