This document describes an implementation of the energy functional minimization proposed by Ballester, Garrido, Lazcano and Caselles for joint optical flow and occlusion estimation. The method build up from the ideas behind the TV-L1 approach introduced by Zach, Pock and Bischof in 2007 but incorporating information that allows to detect occlusions. This information is based on the divergence of the flow and the proposed energy favors the location of occlusions on regions where this divergence is negative. The implemented variational method uses three consecutive frames. The energy functional is composed of regularization terms using the total variation, a data term using the L1 norm, and a term dealing with the occlusions. In the present implementation, we solve the stationary system of partial differential equations arising from the dual minimization problem associated with the TV operator by a variation of the box relaxation numerical scheme proposed by Garamendi, Gaspar, Malpica and Schiavi. This makes the overall algorithm faster than previous implementations based on a gradient descent method.