<p>The Stokes equation is fundamental in fluid mechanics. We used bivariate Bernstein polynomial bases to construct the function space for mixed finite element methods to solve the 2D Stokes equation. Our results show that the numerical accuracy and convergence order using bicubic and lower-order Lagrange interpolation polynomials are comparable to those achieved with Bernstein polynomial bases. However, high-order Lagrange interpolation functions often suffer from the Runge's phenomenon, which limits their effectiveness. By employing high-order Bernstein polynomial bases, we have significantly improved the numerical solutions, effectively mitigating the Runge phenomenon. This approach highlights the advantages of Bernstein polynomial bases in achieving stable and accurate solutions for the 2D Stokes equation.</p>