In this work, we have implemented an orbital-free density functional theory (OF-DFT) solver using the finite element method. In OF-DFT, the total ground state energy is minimized directly with respect to the electron density, rather than via orbitals like in the standard Kohn-Sham approach. For this to be possible, one needs an approximation of a universal density functional of the non-interacting kinetic energy. Presently available approximations allow for computation with very low computational expense, but which gives inaccurate energies. A stable OF-DFT code can be used as a testbed for new kinetic energy functionals and provide the necessary tool for investigating the accuracy of OF-DFT calculations for complex systems. We have implemented Thomas-Fermi theory with and without nuclear cusp condition, as well as additional exchange terms of Dirac and Amaldi. The program uses an extended version of the steepest descent in order to find the minimizing density in the variational principle. Our results include convergence tests for the hydrogen atom, weak bonding in the H2 molecule, and accurate results for the lightest noble gases (He, Ne, Ar). For heavier atoms (Kr, Xe, Rn), the results are less accurate. In addition, we consider hydrogen in the simple cubic structure without the cusp condition, which is a first attempt to use the code for periodic systems. Lastly, we discuss some possible improvements for the iterative process towards the minimizing density, as well as other possible directions for future development.