This work is devoted to the study of a fully discrete scheme for a repulsive chemotaxis with quadratic production model. By following the ideas presented in [Guilen-Gonzalez et al], we introduce an auxiliary variable (the gradient of the chemical concentration), and prove that the corresponding Finite Element (FE) backward Euler scheme is conservative and unconditionally energy-stable. Additionally, we also study some properties like solvability, a priori estimates, convergence towards weak solutions and error estimates. On the other hand, we propose two linear iterative methods to approach the nonlinear scheme: an energy-stable Picard method and Newton's method. We prove solvability and convergence of both methods towards the nonlinear scheme. Finally, we provide some numerical results in agreement with our theoretical analysis with respect to the error estimates.