Summary: In this paper, a difference finite element (DFE) method is presented for the 3D steady Stokes equations. This new method consists of transmitting the finite element solution \(( u_h, p_h)\) of the 3D steady Stokes equations in the direction of \((x, y, z)\) into a series of the finite element solution \(( u_h^k, p_h^k)\) of the 2D steady Stokes equations. Here the 2D steady Stokes equations are solved by the finite element space pair \(( P_1^b, P_1^b, P_1) \times P_1\), where the 2D finite element pair \(( P_1^b, P_1^b) \times P_1\) satisfies the discrete inf-sup condition in a 2D domain \(\omega \). Here we design the weak formulation of the DFE method based on the 3D finite element pair \((( P_1^b, P_1^b, P_1) \times P_1) \times( P_1 \times P_0)\) under the quasi-uniform mesh condition, prove that the 3D finite element pair satisfies the discrete inf-sup condition in a 3D domain \(\Omega\) and provide the existence, uniqueness and stability of the DFE solution \(( u_h, p_h) = ( \sum_{k = 0}^{l_3} u_h^k \phi_k(z), \sum_{k = 1}^{l_3} p_h^k \psi_k(z))\) and deduce the first order convergence of the DFE solution \(( u_h, p_h)\) with respect to the exact solution \((u, p)\) of the 3D steady Stokes equations. Finally, some numerical tests are presented to show the accuracy and efficiency for the proposed method.