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) x P-1, where the 2D finite element pair (P-1(b) , P-1(b)) x 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) x P-1) x ( P-1 x 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)) = ( Sigma(l3)(k=0) u(h)(k) phi(k)(z ), Sigma(l3)(k=1) 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. (C) 2021 IMACS. Published by Elsevier B.V. All rights reserved.