diff --git a/buildconfig/stubs/pygame/draw.pyi b/buildconfig/stubs/pygame/draw.pyi
index 57bae66712..399844b8ba 100644
--- a/buildconfig/stubs/pygame/draw.pyi
+++ b/buildconfig/stubs/pygame/draw.pyi
@@ -89,3 +89,9 @@ def aalines(
     closed: bool,
     points: SequenceLike[Coordinate],
 ) -> Rect: ...
+def bezier(
+    surface: Surface,
+    points: SequenceLike[Coordinate],
+    steps: int,
+    color: ColorLike,
+) -> None: ...
diff --git a/docs/reST/ref/draw.rst b/docs/reST/ref/draw.rst
index 1c71eebe49..3f6afdee2a 100644
--- a/docs/reST/ref/draw.rst
+++ b/docs/reST/ref/draw.rst
@@ -598,6 +598,39 @@ object around the draw calls (see :func:`pygame.Surface.lock` and
 
    .. ## pygame.draw.aalines ##
 
+.. function:: bezier
+
+   | :sl:`draw a Bezier curve`
+   | :sg:`bezier(surface, points, steps, color) -> None`
+
+   Draws a Bézier curve on the given surface.
+
+   :param Surface surface: surface to draw on
+   :param points: a sequence of 3 or more (x, y) coordinates used to form a
+      curve, where each *coordinate* in the sequence must be a
+      tuple/list/:class:`pygame.math.Vector2` of 2 ints/floats (float values
+      will be truncated)
+   :type points: tuple(coordinate) or list(coordinate)
+   :param int steps: number of steps for the interpolation, the minimum is 2
+   :param color: color to draw with, the alpha value is optional if using a
+      tuple ``(RGB[A])``
+   :type color: Color or string (for :doc:`color_list`) or int or tuple(int, int, int, [int])
+
+   :returns: ``None``
+   :rtype: NoneType
+
+   :raises ValueError: if ``steps < 2``
+   :raises ValueError: if ``len(points) < 3`` (must have at least 3 points)
+   :raises IndexError: if ``len(coordinate) < 2`` (each coordinate must have
+      at least 2 items)
+
+   .. note:: This function supports up to around 150-200 points before the algorithm
+             breaks down.
+   
+   .. versionadded:: 2.6.0
+
+   .. ## pygame.draw.bezier ##
+
 .. ## pygame.draw ##
 
 .. figure:: code_examples/draw_module_example.png
diff --git a/src_c/doc/draw_doc.h b/src_c/doc/draw_doc.h
index b47ac37250..4e0f8c69d4 100644
--- a/src_c/doc/draw_doc.h
+++ b/src_c/doc/draw_doc.h
@@ -10,3 +10,4 @@
 #define DOC_DRAW_LINES "lines(surface, color, closed, points) -> Rect\nlines(surface, color, closed, points, width=1) -> Rect\ndraw multiple contiguous straight line segments"
 #define DOC_DRAW_AALINE "aaline(surface, color, start_pos, end_pos) -> Rect\ndraw a straight antialiased line"
 #define DOC_DRAW_AALINES "aalines(surface, color, closed, points) -> Rect\ndraw multiple contiguous straight antialiased line segments"
+#define DOC_DRAW_BEZIER "bezier(surface, points, steps, color) -> None\ndraw a Bezier curve"
diff --git a/src_c/draw.c b/src_c/draw.c
index edd2e55c00..48bd045c7e 100644
--- a/src_c/draw.c
+++ b/src_c/draw.c
@@ -92,6 +92,9 @@ static void
 draw_round_rect(SDL_Surface *surf, int x1, int y1, int x2, int y2, int radius,
                 int width, Uint32 color, int top_left, int top_right,
                 int bottom_left, int bottom_right, int *drawn_area);
+static int
+draw_bezier_color(SDL_Surface *dst, const int *vx, const int *vy, int n, int s,
+                  Uint32 color);
 
 // validation of a draw color
 #define CHECK_LOAD_COLOR(colorobj)                               \
@@ -1105,6 +1108,94 @@ rect(PyObject *self, PyObject *args, PyObject *kwargs)
         return pgRect_New4(rect->x, rect->y, 0, 0);
 }
 
+static PyObject *
+bezier(PyObject *self, PyObject *args, PyObject *kwargs)
+{
+    pgSurfaceObject *surface;
+    PyObject *colorobj, *points, *item;
+    int *vx, *vy, x, y;
+    Py_ssize_t count, i;
+    SDL_Surface *surf = NULL;
+    int ret, steps;
+    int result;
+    Uint32 color;
+
+    static char *keywords[] = {"surface", "points", "steps", "color", NULL};
+
+    // ASSERT_VIDEO_INIT(NULL); To be kept ?
+
+    if (!PyArg_ParseTupleAndKeywords(args, kwargs, "O!OiO", keywords,
+                                     &pgSurface_Type, &surface, &points,
+                                     &steps, &colorobj))
+        return NULL;
+
+    surf = pgSurface_AsSurface(surface);
+    SURF_INIT_CHECK(surf)
+
+    CHECK_LOAD_COLOR(colorobj)
+
+    if (!PySequence_Check(points)) {
+        return RAISE(PyExc_TypeError, "points must be a sequence");
+    }
+
+    count = PySequence_Size(points);
+    if (count < 3) {
+        return RAISE(PyExc_ValueError,
+                     "points must contain more than 2 points");
+    }
+
+    if (steps < 2) {
+        return RAISE(PyExc_ValueError,
+                     "steps parameter must be greater than 1");
+    }
+
+    vx = PyMem_New(int, (size_t)count);
+    vy = PyMem_New(int, (size_t)count);
+    if (!vx || !vy) {
+        if (vx)
+            PyMem_Free(vx);
+        if (vy)
+            PyMem_Free(vy);
+        return RAISE(PyExc_MemoryError, "memory allocation failed");
+    }
+
+    for (i = 0; i < count; i++) {
+        item = PySequence_ITEM(points, i);
+        result = pg_TwoIntsFromObj(item, &x, &y);
+
+        if (!result) {
+            PyMem_Free(vx);
+            PyMem_Free(vy);
+            return RAISE(PyExc_TypeError, "points must be number pairs");
+        }
+
+        Py_DECREF(item);
+
+        vx[i] = x;
+        vy[i] = y;
+    }
+
+    if (!pgSurface_Lock(surface)) {
+        return RAISE(PyExc_RuntimeError, "error locking surface");
+    }
+
+    Py_BEGIN_ALLOW_THREADS;
+    ret = draw_bezier_color(surf, vx, vy, (int)count, steps, color);
+    Py_END_ALLOW_THREADS;
+
+    if (!pgSurface_Unlock(surface)) {
+        return RAISE(PyExc_RuntimeError, "error unlocking surface");
+    }
+
+    PyMem_Free(vx);
+    PyMem_Free(vy);
+
+    if (ret == -1) {
+        return RAISE(pgExc_SDLError, SDL_GetError());
+    }
+    Py_RETURN_NONE;
+}
+
 /* Functions used in drawing algorithms */
 
 static void
@@ -3092,6 +3183,135 @@ draw_round_rect(SDL_Surface *surf, int x1, int y1, int x2, int y2, int radius,
     }
 }
 
+/* ---- Port of Bezier algorithms from SDL_gfxPrimitives.c */
+
+/*!
+\brief Internal function to calculate bezier interpolator of data array with
+ndata values at position 't'.
+
+\param data Array of values.
+\param ndata Size of array.
+\param t Position for which to calculate interpolated value. t should be
+between [0, nstepdata]. \param nstepdata Number of steps for the interpolation
+multiplied by the number of points. \returns Interpolated value at position t,
+value[0] when t<0, value[n-1] when t>=n.
+*/
+static double
+_evaluateBezier(double *data, int ndata, int t, int nstepdata)
+{
+    double mu, result;
+    int n, k, kn, nn, nkn;
+    double blend, muk, munk;
+
+    /* Sanity check bounds */
+    if (t < 0) {
+        return (data[0]);
+    }
+    if (t >= nstepdata) {
+        return (data[ndata - 1]);
+    }
+
+    /* Adjust t to the range 0.0 to 1.0 */
+    mu = t / (double)nstepdata;
+
+    /* Calculate interpolate */
+    n = ndata - 1;
+    result = 0.0;
+    muk = 1;
+    munk = pow(1 - mu, (double)n);
+
+    /* Ensure munk is not 0 which would cause coordinates to be (0, 0) */
+    if (munk <= 0) {
+        return (data[ndata - 1]);
+    }
+
+    for (k = 0; k <= n; k++) {
+        nn = n;
+        kn = k;
+        nkn = n - k;
+        blend = muk * munk;
+        muk *= mu;
+        munk /= (1 - mu);
+        while (nn >= 1) {
+            blend *= nn;
+            nn--;
+            if (kn > 1) {
+                blend /= (double)kn;
+                kn--;
+            }
+            if (nkn > 1) {
+                blend /= (double)nkn;
+                nkn--;
+            }
+        }
+        result += data[k] * blend;
+    }
+
+    return (result);
+}
+
+/*!
+\brief Draw a bezier curve with alpha blending.
+
+\param dst The surface to draw on.
+\param vx Vertex array containing X coordinates of the points of the bezier
+curve. \param vy Vertex array containing Y coordinates of the points of the
+bezier curve. \param n Number of points in the vertex array. Minimum number
+is 3. \param s Number of steps for the interpolation. Minimum number is 2.
+\param color The color value of the bezier curve to draw (0xRRGGBBAA).
+
+\returns Returns 0 on success, -1 on failure.
+*/
+static int
+draw_bezier_color(SDL_Surface *dst, const int *vx, const int *vy, int n, int s,
+                  Uint32 color)
+{
+    int i, steppoints;
+    double *x, *y;
+    int x1, y1, x2, y2;
+    int drawn_area[4] = {INT_MAX, INT_MAX, INT_MIN,
+                         INT_MIN}; /* Used to store bounding box values */
+
+    /*
+     * Variable setup
+     */
+    steppoints = s * n;
+
+    /* Transfer vertices into float arrays */
+    if ((x = (double *)malloc(sizeof(double) * (n + 1))) == NULL) {
+        return (-1);
+    }
+    if ((y = (double *)malloc(sizeof(double) * (n + 1))) == NULL) {
+        free(x);
+        return (-1);
+    }
+    for (i = 0; i < n; i++) {
+        x[i] = (double)vx[i];
+        y[i] = (double)vy[i];
+    }
+    x[n] = (double)vx[0];
+    y[n] = (double)vy[0];
+
+    /*
+     * Draw
+     */
+    x1 = (int)lrint(_evaluateBezier(x, n + 1, 0, steppoints));
+    y1 = (int)lrint(_evaluateBezier(y, n + 1, 0, steppoints));
+    for (i = 1; i <= steppoints; i++) {
+        x2 = (int)_evaluateBezier(x, n, i, steppoints);
+        y2 = (int)_evaluateBezier(y, n, i, steppoints);
+        draw_line(dst, x1, y1, x2, y2, color, drawn_area);
+        x1 = x2;
+        y1 = y2;
+    }
+
+    /* Clean up temporary array */
+    free(x);
+    free(y);
+
+    return 0;
+}
+
 /* List of python functions */
 static PyMethodDef _draw_methods[] = {
     {"aaline", (PyCFunction)aaline, METH_VARARGS | METH_KEYWORDS,
@@ -3111,6 +3331,8 @@ static PyMethodDef _draw_methods[] = {
     {"polygon", (PyCFunction)polygon, METH_VARARGS | METH_KEYWORDS,
      DOC_DRAW_POLYGON},
     {"rect", (PyCFunction)rect, METH_VARARGS | METH_KEYWORDS, DOC_DRAW_RECT},
+    {"bezier", (PyCFunction)bezier, METH_VARARGS | METH_KEYWORDS,
+     DOC_DRAW_BEZIER},
 
     {NULL, NULL, 0, NULL}};
 
diff --git a/test/draw_test.py b/test/draw_test.py
index 3f00a9319e..bad33387ec 100644
--- a/test/draw_test.py
+++ b/test/draw_test.py
@@ -7075,6 +7075,119 @@ class to add any draw.arc specific tests to.
     """
 
 
+### Draw Bezier Testing #######################################################
+
+
+class DrawBezierTest(unittest.TestCase):
+    is_started = False
+
+    foreground_color = (128, 64, 8)
+    background_color = (255, 255, 255)
+
+    def make_palette(base_color):
+        """Return color palette that is various intensities of base_color"""
+        # Need this function for Python 3.x so the base_color
+        # is within the scope of the list comprehension.
+
+        palette = []
+        for i in range(0, 256):
+            r, g, b = base_color
+            if 0 <= i <= 127:
+                # Darken
+                palette.append(((r * i) // 127, (g * i) // 127, (b * i) // 127))
+            else:
+                # Lighten
+                palette.append(
+                    (
+                        r + ((255 - r) * (255 - i)) // 127,
+                        g + ((255 - g) * (255 - i)) // 127,
+                        b + ((255 - b) * (255 - i)) // 127,
+                    )
+                )
+        return palette
+
+    default_palette = make_palette(foreground_color)
+
+    default_size = (100, 100)
+
+    def setUp(self):
+        # This makes sure pygame is always initialized before each test (in
+        # case a test calls pygame.quit()).
+        if not pygame.get_init():
+            pygame.init()
+
+        Surface = pygame.Surface
+        size = self.default_size
+        palette = self.default_palette
+        if not self.is_started:
+            # Create test surfaces
+            self.surfaces = [
+                Surface(size, 0, 8),
+                Surface(size, SRCALPHA, 16),
+                Surface(size, SRCALPHA, 32),
+            ]
+            self.surfaces[0].set_palette(palette)
+            nonpalette_fmts = (
+                # (8, (0xe0, 0x1c, 0x3, 0x0)),
+                (12, (0xF00, 0xF0, 0xF, 0x0)),
+                (15, (0x7C00, 0x3E0, 0x1F, 0x0)),
+                (15, (0x1F, 0x3E0, 0x7C00, 0x0)),
+                (16, (0xF00, 0xF0, 0xF, 0xF000)),
+                (16, (0xF000, 0xF00, 0xF0, 0xF)),
+                (16, (0xF, 0xF0, 0xF00, 0xF000)),
+                (16, (0xF0, 0xF00, 0xF000, 0xF)),
+                (16, (0x7C00, 0x3E0, 0x1F, 0x8000)),
+                (16, (0xF800, 0x7C0, 0x3E, 0x1)),
+                (16, (0x1F, 0x3E0, 0x7C00, 0x8000)),
+                (16, (0x3E, 0x7C0, 0xF800, 0x1)),
+                (16, (0xF800, 0x7E0, 0x1F, 0x0)),
+                (16, (0x1F, 0x7E0, 0xF800, 0x0)),
+                (24, (0xFF, 0xFF00, 0xFF0000, 0x0)),
+                (24, (0xFF0000, 0xFF00, 0xFF, 0x0)),
+                (32, (0xFF0000, 0xFF00, 0xFF, 0x0)),
+                (32, (0xFF000000, 0xFF0000, 0xFF00, 0x0)),
+                (32, (0xFF, 0xFF00, 0xFF0000, 0x0)),
+                (32, (0xFF00, 0xFF0000, 0xFF000000, 0x0)),
+                (32, (0xFF0000, 0xFF00, 0xFF, 0xFF000000)),
+                (32, (0xFF000000, 0xFF0000, 0xFF00, 0xFF)),
+                (32, (0xFF, 0xFF00, 0xFF0000, 0xFF000000)),
+                (32, (0xFF00, 0xFF0000, 0xFF000000, 0xFF)),
+            )
+            for bitsize, masks in nonpalette_fmts:
+                self.surfaces.append(Surface(size, 0, bitsize, masks))
+        for surf in self.surfaces:
+            surf.fill(self.background_color)
+
+    def check_at(self, surf, posn, color):
+        sc = surf.get_at(posn)
+        depth = surf.get_bitsize()
+        flags = surf.get_flags()
+        masks = surf.get_masks()
+        fail_msg = f"{sc} != {color} at {posn}, bitsize: {depth}, flags: {flags}, masks: {masks}"
+        self.assertEqual(sc, color, fail_msg)
+
+    def test_bezier(self):
+        """bezier(surface, points, steps, color): return None"""
+        fg = self.foreground_color
+        bg = self.background_color
+        points = [(10, 50), (25, 15), (60, 80), (92, 30)]
+        fg_test_points = [points[0], points[3]]
+        bg_test_points = [
+            (points[0][0] - 1, points[0][1]),
+            (points[3][0] + 1, points[3][1]),
+            (points[1][0], points[1][1] + 3),
+            (points[2][0], points[2][1] - 3),
+        ]
+        for surf in self.surfaces:
+            fg_adjusted = surf.unmap_rgb(surf.map_rgb(fg))
+            bg_adjusted = surf.unmap_rgb(surf.map_rgb(bg))
+            pygame.draw.bezier(surf, points, 30, fg)
+            for posn in fg_test_points:
+                self.check_at(surf, posn, fg_adjusted)
+            for posn in bg_test_points:
+                self.check_at(surf, posn, bg_adjusted)
+
+
 ### Draw Module Testing #######################################################