We develop the first algorithm able to jointly compute the maximum a posteriori estimate of the Cosmic Microwave Background (CMB) temperature and polarization fields, the gravitational potential by which they are lensed, and the tensor-to-scalar ratio, r . This is an important step towards sampling from the joint posterior probability function of these quantities, which, assuming Gaussianity of the CMB fields and lensing potential, contains all available cosmological information and would yield theoretically optimal constraints. Attaining such optimal constraints will be crucial for next-generation CMB surveys like CMB-S4, where limits on r could be improved by factors of a few over currently used suboptimal quadratic estimators. The maximization procedure described here depends on a newly developed lensing algorithm, which we term LenseFlow, and which lenses a map by solving a system of ordinary differential equations. This description has conceptual advantages, such as allowing us to give a simple nonperturbative proof that the determinant of LenseFlow on pixelized maps is equal to unity, which is crucial for our purposes and unique to LenseFlow as compared to other lensing algorithms we have tested. It also has other useful properties such as that it can be trivially inverted (i.e., delensing) for the same computational cost as the forward operation, and can be used to compute lensing adjoint, Jacobian, and Hessian operators. We test and validate the maximization procedure on flat-sky simulations covering up to 600 deg2 with nonwhite noise and masking.