En este trabajo, estudiamos un modelo de remodelación ósea utilizado para reproducir el fenómeno de la osteointegración alrededor de implantes endoóseos. El problema biológico está escrito en términos de las densidades de plaquetas, células osteogénicas y osteoblastos y las concentraciones de dos factores de crecimiento. Su formulación variacional conduce a un sistema fuertemente acoplado de ecuaciones variacionales parabólicas no lineales. Se establece un resultado de existencia y unicidad de esta forma variacional. Luego, se introduce una aproximación completamente discreta del problema utilizando el método de elementos finitos y un esquema de Euler semi-implícito. Se obtienen estimaciones de error a priori y se deriva la convergencia lineal del algoritmo bajo algunas condiciones de regularidad adecuadas y se prueba con un ejemplo numérico. Finalmente, se presentan resultados numéricos de una y dos dimensiones para demostrar la precisión del algoritmo y el comportamiento de la solución.