This study addresses how to implement the Galerkin finite element and least square finite element methods using auxiliary equations to solve the partial differential equation numerically, which models the convection-diffusion-reaction, set on a steady 3D domain. In the spatial discretization, hexahedral elements with eight (linear element) and 27 (quadratic element) nodes were used, in which Lagrange interpolation functions were adopted in local coordinates. Turning all the formulation of the problem of global coordinates into local coordinates, the Gauss-Legendre quadrature method was used to integrate coefficients of the element matrices numerically. In addition to the formulation by the two methods, a computer code was implemented to simulate the phenomenon proposed. By using analytical solutions, sundry numerical error analysis was performed from L-2 norm (domain-average error) and L-infinity norm (domain-top error), thus validating the numerical results. A real case is proposed and assessed.