The flexo-photovoltaic effect has emerged as a promising alternative to classical photovoltaic technologies, taking advantage of flexoelectric-induced electric fields to generate charge separation. While early studies in flexotronics focused on carrier redistribution without charge transport or relied on linearized transport equations, most existing models are restricted to simplified geometries using small perturbation theories. Fully coupled, nonlinear continuum models incorporating charge transport remain largely unexplored, particularly within a Finite Element (FE) framework. This work addresses this gap by presenting the FE solution of the fourth-order, nonlinear system of partial differential equations modeling flexo-photovoltaics at infinitesimal deformations, without resorting to linearization assumptions, by means of a C0 interior penalty method. A three-dimensional flexotronic plate under bending is simulated, and an illuminated p-n junction device is generalized from conventional photovoltaics to flexo-photovoltaics, demonstrating the ability to reproduce standard photovoltaic device behavior through flexoelectricity. Additionally, a continuum model for flexo-photovoltaics at finite deformations is proposed, incorporating drift effects arising from volume changes in charge transport. The results contribute to the numerical and theoretical foundations of flexo-photovoltaics and highlight the potential of flexoelectricity for next-generation energy conversion devices.