In this work, we study a bone remodeling model used to reproduce the phenomenon of osseointegration around endosseous implants. The biological problem is written in terms of the densities of platelets, osteogenic cells, and osteoblasts and the concentrations of two growth factors. Its variational formulation leads to a strongly coupled nonlinear system of parabolic variational equations. An existence and uniqueness result of this variational form is stated. Then, a fully discrete approximation of the problem is introduced by using the finite element method and a semi-implicit Euler scheme. A priori error estimates are obtained, and the linear convergence of the algorithm is derived under some suitable regularity conditions and tested with a numerical example. Finally, one- and two-dimensional numerical results are presented to demonstrate the accuracy of the algorithm and the behavior of the solution.