A modified, hybrid nodal-integral/finite-element method (NI-FEM) is developed to solve the three-dimensional (3D), steady-state, convection-diffusion problems in arbitrary geometries. The hybrid NI-FEM takes advantage of the high efficiency of the conventional nodal-integral method (NIM) and the flexible mesh generation of the finite element method (FEM), which is applicable to arbitrary geometries. In this method, the computational domain is discretized into, for 3D problems, four-node tetrahedral elements and eight-node cuboid elements. The cuboid elements are used to discretize the interior region and regions adjacent to boundaries that are parallel to planes formed by two of the axes, while the tetrahedral elements are used to discretize the remaining irregular regions. The conventional NIM is used to develop difference equations for the transverse-averaged variables on the interface between two adjacent cuboid elements. The FEM is applied to develop algebraic equations for the node temperatures of tetrahedral elements. On the interface between two different kinds of elements, the transverse averaged variable for the cuboid element is obtained by averaging the node values of its adjacent tetrahedral elements, while the heat flux for the tetrahedral element is calculated using the corresponding transverse-averaged variables of its adjacent cuboid element. The hybrid NI-FEM is developed to be solved using a matrix formulation for the entire domain rather than an iterative procedure, detailed derivations of which for the 3D case are presented in this paper. Using the NI-FEM developed here, 2D and 3D convection-diffusion test problems are solved, and the numerical results are compared to the exact (manufactured) solutions to determine the order, accuracy and efficiency of the method. Numerical scheme is found to be of, as expected, second order. (C) 2018 Elsevier Ltd. All rights reserved.