diff --git a/docs/source/reference.rst b/docs/source/reference.rst
index 5352f83223..771c222922 100644
--- a/docs/source/reference.rst
+++ b/docs/source/reference.rst
@@ -67,6 +67,7 @@ Mobjects
    manim.mobject.matrix
    manim.mobject.mobject
    manim.mobject.table
+   manim.mobject.three_d.implicit_surface
    manim.mobject.three_d.polyhedra
    manim.mobject.three_d.three_d_utils
    manim.mobject.three_d.three_dimensions
diff --git a/example_scenes/implicit_surface_examples.py b/example_scenes/implicit_surface_examples.py
new file mode 100644
index 0000000000..aa43c0c3e6
--- /dev/null
+++ b/example_scenes/implicit_surface_examples.py
@@ -0,0 +1,185 @@
+"""Example scenes demonstrating ImplicitSurface usage."""
+
+from __future__ import annotations
+
+import numpy as np
+
+from manim import *
+
+
+class ImplicitSphereExample(ThreeDScene):
+    """Demonstrate a unit sphere defined implicitly by x^2 + y^2 + z^2 = 1."""
+
+    def construct(self):
+        axes = ThreeDAxes(
+            x_range=(-1.5, 1.5, 1),
+            y_range=(-1.5, 1.5, 1),
+            z_range=(-1.5, 1.5, 1),
+        )
+
+        def sphere_func(x, y, z):
+            return x**2 + y**2 + z**2 - 1.0
+
+        surface = ImplicitSurface(
+            sphere_func,
+            x_range=(-1.3, 1.3),
+            y_range=(-1.3, 1.3),
+            z_range=(-1.3, 1.3),
+            resolution=40,
+            fill_color=BLUE,
+            fill_opacity=0.8,
+        )
+
+        self.set_camera_orientation(phi=70 * DEGREES, theta=-45 * DEGREES)
+        self.add(axes, surface)
+        self.wait()
+
+
+class ImplicitTorusExample(ThreeDScene):
+    """Demonstrate a torus defined implicitly."""
+
+    def construct(self):
+        R, r = 2.0, 0.5  # Major and minor radii
+
+        def torus_func(x, y, z):
+            return (np.sqrt(x**2 + y**2) - R) ** 2 + z**2 - r**2
+
+        surface = ImplicitSurface(
+            torus_func,
+            x_range=(-3, 3),
+            y_range=(-3, 3),
+            z_range=(-1, 1),
+            resolution=50,
+            fill_color=GREEN,
+            fill_opacity=0.9,
+        )
+
+        self.set_camera_orientation(phi=60 * DEGREES, theta=-45 * DEGREES)
+        self.add(surface)
+        self.wait()
+
+
+class ImplicitGyroidExample(ThreeDScene):
+    """Demonstrate a gyroid minimal surface with gradient coloring."""
+
+    def construct(self):
+        def gyroid_func(x, y, z):
+            return np.sin(x) * np.cos(y) + np.sin(y) * np.cos(z) + np.sin(z) * np.cos(x)
+
+        # Define salmon to pink gradient colors
+        salmon = ManimColor("#FA8072")
+        hot_pink = ManimColor("#FF69B4")
+
+        surface = ImplicitSurface(
+            gyroid_func,
+            x_range=(-PI, PI),
+            y_range=(-PI, PI),
+            z_range=(-PI, PI),
+            resolution=50,
+            fill_opacity=0.9,
+        )
+
+        # Apply gradient coloring based on z-position
+        for tri in surface.submobjects:
+            # Get the center z-coordinate of each triangle
+            z_val = tri.get_center()[2]
+            # Normalize z to [0, 1] range
+            t = (z_val + PI) / (2 * PI)
+            t = max(0, min(1, t))  # Clamp to [0, 1]
+            # Interpolate from salmon to pink
+            tri.set_fill(
+                color=interpolate_color(salmon, hot_pink, t),
+                opacity=0.9,
+            )
+
+        self.set_camera_orientation(phi=70 * DEGREES, theta=-30 * DEGREES)
+        self.add(surface)
+        self.wait()
+
+
+class ImplicitBlobbySurfaceExample(ThreeDScene):
+    """Demonstrate a 'blobby' surface made from sum of Gaussians."""
+
+    def construct(self):
+        centers = [
+            (-0.5, 0.0, 0.0),
+            (0.5, 0.0, 0.0),
+            (0.0, 0.8, 0.0),
+        ]
+
+        def blobby_func(x, y, z):
+            # Sum of Gaussians centered at different points
+            field = 0.0
+            for cx, cy, cz in centers:
+                r2 = (x - cx) ** 2 + (y - cy) ** 2 + (z - cz) ** 2
+                field += np.exp(-4 * r2)
+            return field - 0.4  # 0-level of this isosurface
+
+        surface = ImplicitSurface(
+            blobby_func,
+            x_range=(-2, 2),
+            y_range=(-2, 2),
+            z_range=(-2, 2),
+            resolution=40,
+            fill_color=ORANGE,
+            fill_opacity=0.8,
+        )
+
+        axes = ThreeDAxes(
+            x_range=(-2, 2, 1),
+            y_range=(-2, 2, 1),
+            z_range=(-2, 2, 1),
+        )
+
+        self.set_camera_orientation(phi=70 * DEGREES, theta=-60 * DEGREES)
+        self.add(axes, surface)
+        self.wait()
+
+
+class ImplicitCylinderExample(ThreeDScene):
+    """Demonstrate a cylinder defined implicitly by x^2 + y^2 = 1."""
+
+    def construct(self):
+        def cylinder_func(x, y, z):
+            return x**2 + y**2 - 1.0
+
+        surface = ImplicitSurface(
+            cylinder_func,
+            x_range=(-1.5, 1.5),
+            y_range=(-1.5, 1.5),
+            z_range=(-1.5, 1.5),
+            resolution=40,
+            fill_color=YELLOW,
+            fill_opacity=0.7,
+        )
+
+        self.set_camera_orientation(phi=70 * DEGREES, theta=-45 * DEGREES)
+        self.add(surface)
+        self.wait()
+
+
+class ImplicitHeartExample(ThreeDScene):
+    """Demonstrate a 3D heart shape defined implicitly."""
+
+    def construct(self):
+        def heart_func(x, y, z):
+            # Heart surface equation
+            return (
+                (x**2 + (9 / 4) * y**2 + z**2 - 1) ** 3
+                - x**2 * z**3
+                - (9 / 80) * y**2 * z**3
+            )
+
+        surface = ImplicitSurface(
+            heart_func,
+            x_range=(-1.5, 1.5),
+            y_range=(-1.5, 1.5),
+            z_range=(-1.5, 1.5),
+            resolution=50,
+            fill_color=RED,
+            fill_opacity=0.9,
+        )
+
+        self.set_camera_orientation(phi=70 * DEGREES, theta=-45 * DEGREES)
+        self.add(surface)
+        self.wait()
diff --git a/manim/__init__.py b/manim/__init__.py
index 0605d4a3ae..e5a4d071aa 100644
--- a/manim/__init__.py
+++ b/manim/__init__.py
@@ -73,6 +73,7 @@
 from .mobject.text.numbers import *
 from .mobject.text.tex_mobject import *
 from .mobject.text.text_mobject import *
+from .mobject.three_d.implicit_surface import *
 from .mobject.three_d.polyhedra import *
 from .mobject.three_d.three_d_utils import *
 from .mobject.three_d.three_dimensions import *
diff --git a/manim/mobject/opengl/opengl_compatibility.py b/manim/mobject/opengl/opengl_compatibility.py
index 94991d544e..e5ddae6d8e 100644
--- a/manim/mobject/opengl/opengl_compatibility.py
+++ b/manim/mobject/opengl/opengl_compatibility.py
@@ -3,6 +3,7 @@
 from abc import ABCMeta
 
 from manim import config
+from manim.mobject.opengl.opengl_implicit_surface import OpenGLImplicitSurface
 from manim.mobject.opengl.opengl_mobject import OpenGLMobject
 from manim.mobject.opengl.opengl_point_cloud_mobject import OpenGLPMobject
 from manim.mobject.opengl.opengl_three_dimensions import OpenGLSurface
@@ -32,6 +33,7 @@ def __new__(mcls, name, bases, namespace):
                 "Mobject1D": OpenGLPMobject,
                 "Mobject2D": OpenGLPMobject,
                 "Surface": OpenGLSurface,
+                "ImplicitSurface": OpenGLImplicitSurface,
             }
 
             bases = tuple(
diff --git a/manim/mobject/opengl/opengl_implicit_surface.py b/manim/mobject/opengl/opengl_implicit_surface.py
new file mode 100644
index 0000000000..f965a6e5d3
--- /dev/null
+++ b/manim/mobject/opengl/opengl_implicit_surface.py
@@ -0,0 +1,227 @@
+"""OpenGL implementation of implicit surfaces."""
+
+from __future__ import annotations
+
+from collections.abc import Callable, Sequence
+from typing import Any
+
+import mcubes
+import moderngl
+import numpy as np
+
+from manim.mobject.opengl.opengl_mobject import OpenGLMobject
+from manim.utils.color import GREY, ParsableManimColor
+
+__all__ = ["OpenGLImplicitSurface"]
+
+
+class OpenGLImplicitSurface(OpenGLMobject):
+    """OpenGL implementation of a 3D implicit surface.
+
+    This class creates a triangular mesh representation of an implicit surface
+    using the marching cubes algorithm, optimized for OpenGL rendering.
+
+    Parameters
+    ----------
+    func
+        A callable that takes three arguments (x, y, z) and returns a scalar value.
+        The function should be NumPy-aware for efficient evaluation over grids.
+    x_range
+        The range of x values as (x_min, x_max). Defaults to (-1.0, 1.0).
+    y_range
+        The range of y values as (y_min, y_max). Defaults to (-1.0, 1.0).
+    z_range
+        The range of z values as (z_min, z_max). Defaults to (-1.0, 1.0).
+    resolution
+        The number of sample points along each axis. Can be an integer
+        (same resolution for all axes) or a 3-tuple (nx, ny, nz).
+        Defaults to 32.
+    level
+        The isosurface level to extract. Defaults to 0.0.
+    color
+        The color of the surface. Defaults to GREY.
+    opacity
+        The opacity of the surface. Defaults to 1.0.
+    gloss
+        The glossiness of the surface. Defaults to 0.3.
+    shadow
+        The shadow intensity. Defaults to 0.4.
+    """
+
+    shader_dtype = [
+        ("point", np.float32, (3,)),
+        ("normal", np.float32, (3,)),
+        ("color", np.float32, (4,)),
+    ]
+    shader_folder = "surface"
+
+    def __init__(
+        self,
+        func: Callable[[float, float, float], float] | None = None,
+        x_range: Sequence[float] = (-1.0, 1.0),
+        y_range: Sequence[float] = (-1.0, 1.0),
+        z_range: Sequence[float] = (-1.0, 1.0),
+        resolution: int | Sequence[int] = 32,
+        level: float = 0.0,
+        color: ParsableManimColor = GREY,
+        opacity: float = 1.0,
+        gloss: float = 0.3,
+        shadow: float = 0.4,
+        render_primitive: int = moderngl.TRIANGLES,
+        depth_test: bool = True,
+        **kwargs: Any,
+    ) -> None:
+        self.func = func
+        self.x_range = tuple(x_range)
+        self.y_range = tuple(y_range)
+        self.z_range = tuple(z_range)
+        self.level = level
+        self.resolution = self._normalize_resolution(resolution)
+
+        # Store mesh data
+        self._mesh_vertices: np.ndarray | None = None
+        self._mesh_triangles: np.ndarray | None = None
+        self._triangle_indices: np.ndarray | None = None
+
+        super().__init__(
+            color=color,
+            opacity=opacity,
+            gloss=gloss,
+            shadow=shadow,
+            render_primitive=render_primitive,
+            depth_test=depth_test,
+            **kwargs,
+        )
+
+    @staticmethod
+    def _normalize_resolution(resolution: int | Sequence[int]) -> tuple[int, int, int]:
+        """Convert resolution to a 3-tuple of integers."""
+        if isinstance(resolution, int):
+            if resolution < 2:
+                raise ValueError("resolution must be >= 2")
+            return (resolution, resolution, resolution)
+
+        res = tuple(resolution)
+        if len(res) == 1:
+            if res[0] < 2:
+                raise ValueError("resolution must be >= 2")
+            return (res[0], res[0], res[0])
+        if len(res) != 3:
+            raise ValueError("resolution must be an int or a 3-tuple")
+        if any(r < 2 for r in res):
+            raise ValueError("each resolution component must be >= 2")
+        return (res[0], res[1], res[2])
+
+    def init_points(self) -> None:
+        """Initialize the surface mesh using marching cubes."""
+        if self.func is None:
+            self.set_points(np.zeros((0, 3)))
+            self._triangle_indices = np.zeros(0, dtype=int)
+            return
+
+        nx, ny, nz = self.resolution
+        x_min, x_max = self.x_range
+        y_min, y_max = self.y_range
+        z_min, z_max = self.z_range
+
+        # Create grid
+        xs = np.linspace(x_min, x_max, nx)
+        ys = np.linspace(y_min, y_max, ny)
+        zs = np.linspace(z_min, z_max, nz)
+        X, Y, Z = np.meshgrid(xs, ys, zs, indexing="ij")
+
+        # Evaluate function
+        values = np.asarray(self.func(X, Y, Z), dtype=np.float64)
+
+        # Check if there's a surface
+        if np.all(values > self.level) or np.all(values < self.level):
+            self.set_points(np.zeros((0, 3)))
+            self._triangle_indices = np.zeros(0, dtype=int)
+            return
+
+        try:
+            vertices, triangles = mcubes.marching_cubes(values, self.level)
+        except Exception:
+            self.set_points(np.zeros((0, 3)))
+            self._triangle_indices = np.zeros(0, dtype=int)
+            return
+
+        if len(vertices) == 0:
+            self.set_points(np.zeros((0, 3)))
+            self._triangle_indices = np.zeros(0, dtype=int)
+            return
+
+        # Map to world coordinates
+        vertices[:, 0] = x_min + vertices[:, 0] * (x_max - x_min) / (nx - 1)
+        vertices[:, 1] = y_min + vertices[:, 1] * (y_max - y_min) / (ny - 1)
+        vertices[:, 2] = z_min + vertices[:, 2] * (z_max - z_min) / (nz - 1)
+
+        self._mesh_vertices = vertices
+        self._mesh_triangles = triangles
+
+        # Flatten triangles for rendering
+        all_verts = vertices[triangles.flatten()]
+        self._triangle_indices = np.arange(len(all_verts), dtype=int)
+
+        self.set_points(all_verts)
+
+    def get_triangle_indices(self) -> np.ndarray:
+        """Return indices for triangle rendering."""
+        if self._triangle_indices is None:
+            return np.zeros(0, dtype=int)
+        return self._triangle_indices
+
+    def get_shader_vert_indices(self) -> Sequence[int] | None:
+        """Return vertex indices for shader."""
+        indices = self.get_triangle_indices()
+        return indices.tolist() if len(indices) > 0 else None
+
+    def _compute_normals(self) -> np.ndarray:
+        """Compute normals for each vertex based on face normals."""
+        points = self.points
+        if len(points) == 0:
+            return np.zeros((0, 3))
+
+        normals = np.zeros_like(points)
+
+        # Process triangles (every 3 points is a triangle)
+        for i in range(0, len(points), 3):
+            if i + 2 >= len(points):
+                break
+            p0, p1, p2 = points[i], points[i + 1], points[i + 2]
+            # Compute face normal
+            v1 = p1 - p0
+            v2 = p2 - p0
+            normal = np.cross(v1, v2)
+            norm = np.linalg.norm(normal)
+            if norm > 0:
+                normal = normal / norm
+            # Assign same normal to all 3 vertices
+            normals[i] = normal
+            normals[i + 1] = normal
+            normals[i + 2] = normal
+
+        return normals
+
+    def get_shader_data(self) -> np.ndarray:
+        """Generate shader data for rendering."""
+        points = self.points
+        if len(points) == 0:
+            return np.zeros(0, dtype=self.shader_dtype)
+
+        shader_data = np.zeros(len(points), dtype=self.shader_dtype)
+        shader_data["point"] = points
+
+        # Compute and add normals
+        normals = self._compute_normals()
+        shader_data["normal"] = normals
+
+        # Add colors
+        self.fill_in_shader_color_info(shader_data)
+
+        return shader_data
+
+    def fill_in_shader_color_info(self, shader_data: np.ndarray) -> np.ndarray:
+        """Fill color information in shader data."""
+        self.read_data_to_shader(shader_data, "color", "rgbas")
+        return shader_data
diff --git a/manim/mobject/three_d/__init__.py b/manim/mobject/three_d/__init__.py
index 98d295a24e..1abe7b1877 100644
--- a/manim/mobject/three_d/__init__.py
+++ b/manim/mobject/three_d/__init__.py
@@ -6,6 +6,7 @@
 .. autosummary::
     :toctree: ../reference
 
+    ~implicit_surface
     ~polyhedra
     ~three_d_utils
     ~three_dimensions
diff --git a/manim/mobject/three_d/implicit_surface.py b/manim/mobject/three_d/implicit_surface.py
new file mode 100644
index 0000000000..2dec1e7760
--- /dev/null
+++ b/manim/mobject/three_d/implicit_surface.py
@@ -0,0 +1,315 @@
+"""Three-dimensional implicit surfaces."""
+
+from __future__ import annotations
+
+__all__ = ["ImplicitSurface"]
+
+from collections.abc import Callable, Sequence
+from typing import TYPE_CHECKING, Any
+
+import mcubes
+import numpy as np
+
+from manim.mobject.opengl.opengl_compatibility import ConvertToOpenGL
+from manim.mobject.three_d.three_dimensions import ThreeDVMobject
+from manim.mobject.types.vectorized_mobject import VGroup
+from manim.utils.color import BLUE, ManimColor, ParsableManimColor
+
+if TYPE_CHECKING:
+    pass
+
+
+class ImplicitSurface(VGroup, metaclass=ConvertToOpenGL):
+    """A 3D implicit surface defined by f(x, y, z) = 0.
+
+    This class creates a triangular mesh representation of an implicit surface
+    using the marching cubes algorithm. The surface is defined by a scalar
+    function f(x, y, z), and the class extracts the isosurface where f = level.
+
+    Parameters
+    ----------
+    func
+        A callable that takes three arguments (x, y, z) and returns a scalar value.
+        The function should be NumPy-aware for efficient evaluation over grids.
+        The surface is extracted where func(x, y, z) = level.
+    x_range
+        The range of x values as (x_min, x_max). Defaults to (-1.0, 1.0).
+    y_range
+        The range of y values as (y_min, y_max). Defaults to (-1.0, 1.0).
+    z_range
+        The range of z values as (z_min, z_max). Defaults to (-1.0, 1.0).
+    resolution
+        The number of sample points along each axis. Can be an integer
+        (same resolution for all axes) or a 3-tuple (nx, ny, nz).
+        Higher values give smoother surfaces but take longer to compute.
+        Defaults to 32.
+    level
+        The isosurface level to extract. The surface is where
+        func(x, y, z) = level. Defaults to 0.0.
+    fill_color
+        The fill color of the surface faces. Defaults to BLUE.
+    fill_opacity
+        The opacity of the surface faces, from 0 (transparent) to 1 (opaque).
+        Defaults to 1.0.
+    stroke_color
+        The stroke color of the triangle edges. If None, no stroke is applied.
+    stroke_width
+        The width of the triangle edge strokes. Defaults to 0.0 (no stroke).
+    faces_config
+        Additional keyword arguments passed to each triangle face.
+
+    Examples
+    --------
+    .. manim:: ImplicitSphereExample
+        :save_last_frame:
+
+        class ImplicitSphereExample(ThreeDScene):
+            def construct(self):
+                axes = ThreeDAxes(
+                    x_range=(-1.5, 1.5, 1),
+                    y_range=(-1.5, 1.5, 1),
+                    z_range=(-1.5, 1.5, 1),
+                )
+
+                # Unit sphere: x^2 + y^2 + z^2 = 1
+                def sphere_func(x, y, z):
+                    return x**2 + y**2 + z**2 - 1.0
+
+                surface = ImplicitSurface(
+                    sphere_func,
+                    x_range=(-1.3, 1.3),
+                    y_range=(-1.3, 1.3),
+                    z_range=(-1.3, 1.3),
+                    resolution=40,
+                    fill_color=BLUE,
+                    fill_opacity=0.8,
+                )
+
+                self.set_camera_orientation(phi=70 * DEGREES, theta=-45 * DEGREES)
+                self.add(axes, surface)
+
+    .. manim:: ImplicitTorusExample
+        :save_last_frame:
+
+        class ImplicitTorusExample(ThreeDScene):
+            def construct(self):
+                # Torus with major radius 2, minor radius 0.5
+                R, r = 2.0, 0.5
+
+                def torus_func(x, y, z):
+                    return (np.sqrt(x**2 + y**2) - R)**2 + z**2 - r**2
+
+                surface = ImplicitSurface(
+                    torus_func,
+                    x_range=(-3, 3),
+                    y_range=(-3, 3),
+                    z_range=(-1, 1),
+                    resolution=50,
+                    fill_color=GREEN,
+                    fill_opacity=0.9,
+                )
+
+                self.set_camera_orientation(phi=60 * DEGREES, theta=-45 * DEGREES)
+                self.add(surface)
+
+    Notes
+    -----
+    The function should be NumPy-aware and support broadcasting. If your function
+    only accepts scalar inputs, you can wrap it with ``np.vectorize``::
+
+        def scalar_f(x, y, z):
+            # some computation that only works with scalars
+            return result
+
+
+        # Wrap for use with ImplicitSurface
+        surface = ImplicitSurface(np.vectorize(scalar_f), ...)
+
+    However, vectorized functions are slower than native NumPy implementations.
+
+    See Also
+    --------
+    :class:`.ImplicitFunction` : 2D implicit curves defined by f(x, y) = 0.
+    :class:`.Surface` : Parametric surfaces defined by (u, v) -> (x, y, z).
+    """
+
+    def __init__(
+        self,
+        func: Callable[[float, float, float], float],
+        x_range: Sequence[float] = (-1.0, 1.0),
+        y_range: Sequence[float] = (-1.0, 1.0),
+        z_range: Sequence[float] = (-1.0, 1.0),
+        resolution: int | Sequence[int] = 32,
+        level: float = 0.0,
+        fill_color: ParsableManimColor = BLUE,
+        fill_opacity: float = 1.0,
+        stroke_color: ParsableManimColor | None = None,
+        stroke_width: float = 0.0,
+        faces_config: dict[str, Any] | None = None,
+        **kwargs: Any,
+    ) -> None:
+        super().__init__(**kwargs)
+
+        self.func = func
+        self.x_range = tuple(x_range)
+        self.y_range = tuple(y_range)
+        self.z_range = tuple(z_range)
+        self.level = level
+        self.resolution = self._normalize_resolution(resolution)
+
+        # Build face style configuration
+        self._fill_color = ManimColor(fill_color)
+        self._fill_opacity = fill_opacity
+        self._stroke_color = ManimColor(stroke_color) if stroke_color else None
+        self._stroke_width = stroke_width
+        self._faces_config = faces_config or {}
+
+        self._build_surface()
+
+    @staticmethod
+    def _normalize_resolution(resolution: int | Sequence[int]) -> tuple[int, int, int]:
+        """Convert resolution to a 3-tuple of integers."""
+        if isinstance(resolution, int):
+            if resolution < 2:
+                raise ValueError("resolution must be >= 2")
+            return (resolution, resolution, resolution)
+
+        res = tuple(resolution)
+        if len(res) == 1:
+            if res[0] < 2:
+                raise ValueError("resolution must be >= 2")
+            return (res[0], res[0], res[0])
+        if len(res) != 3:
+            raise ValueError("resolution must be an int or a 3-tuple")
+        if any(r < 2 for r in res):
+            raise ValueError("each resolution component must be >= 2")
+        return (res[0], res[1], res[2])
+
+    def _build_surface(self) -> None:
+        """Build the triangular mesh from the implicit function."""
+        nx, ny, nz = self.resolution
+
+        x_min, x_max = self.x_range
+        y_min, y_max = self.y_range
+        z_min, z_max = self.z_range
+
+        # Create grid coordinates
+        xs = np.linspace(x_min, x_max, nx)
+        ys = np.linspace(y_min, y_max, ny)
+        zs = np.linspace(z_min, z_max, nz)
+
+        # Create meshgrid with indexing="ij" so axis 0 -> x, axis 1 -> y, axis 2 -> z
+        X, Y, Z = np.meshgrid(xs, ys, zs, indexing="ij")
+
+        # Evaluate the function on the grid
+        values = np.asarray(self.func(X, Y, Z), dtype=np.float64)
+        if values.shape != (nx, ny, nz):
+            raise ValueError(
+                f"func must return array of shape {(nx, ny, nz)}, "
+                f"got {values.shape} instead. Ensure your function is NumPy-aware "
+                "and supports broadcasting."
+            )
+
+        # Check if there's actually a surface to extract
+        if np.all(values > self.level) or np.all(values < self.level):
+            # No surface crosses through the domain
+            return
+
+        # Run marching cubes algorithm using PyMCubes
+        try:
+            vertices, triangles = mcubes.marching_cubes(values, self.level)
+        except Exception:
+            # marching_cubes can fail if no surface is found
+            return
+
+        if len(vertices) == 0 or len(triangles) == 0:
+            return
+
+        # Map grid indices in vertices to world coordinates
+        # vertices[:, 0] is in [0, nx-1], etc.
+        vertices[:, 0] = x_min + vertices[:, 0] * (x_max - x_min) / (nx - 1)
+        vertices[:, 1] = y_min + vertices[:, 1] * (y_max - y_min) / (ny - 1)
+        vertices[:, 2] = z_min + vertices[:, 2] * (z_max - z_min) / (nz - 1)
+
+        # Store mesh data
+        self._mesh_vertices = vertices
+        self._mesh_triangles = triangles
+
+        # Build triangle faces using ThreeDVMobject
+        for tri in triangles:
+            p1, p2, p3 = vertices[tri]
+            # Create a closed triangle path: p1 -> p2 -> p3 -> p1
+            face = ThreeDVMobject()
+            face.set_points_as_corners([p1, p2, p3, p1])
+            face.set_fill(
+                color=self._fill_color,
+                opacity=self._fill_opacity,
+            )
+            face.set_stroke(
+                color=self._stroke_color,
+                width=self._stroke_width,
+            )
+            self.add(face)
+
+    def get_vertices(self) -> np.ndarray:
+        """Return all unique vertices of the surface mesh.
+
+        Returns
+        -------
+        np.ndarray
+            An array of shape (N, 3) containing the 3D coordinates
+            of all unique vertices in the mesh.
+        """
+        if hasattr(self, "_mesh_vertices") and self._mesh_vertices is not None:
+            return self._mesh_vertices.copy()
+
+        return np.array([]).reshape(0, 3)
+
+    def verify_surface(
+        self,
+        n_samples: int = 500,
+    ) -> dict[str, float]:
+        """Verify that vertices lie on the implicit surface.
+
+        This method samples vertices from the surface and checks how close
+        they are to satisfying f(x, y, z) = level.
+
+        Parameters
+        ----------
+        n_samples
+            Number of vertices to sample for verification. Defaults to 500.
+
+        Returns
+        -------
+        dict
+            A dictionary with 'mean_error' and 'max_error' keys containing
+            the mean and maximum absolute deviation from the level value.
+
+        Examples
+        --------
+        >>> def sphere(x, y, z):
+        ...     return x**2 + y**2 + z**2 - 1
+        >>> surface = ImplicitSurface(sphere, resolution=20)
+        >>> errors = surface.verify_surface()
+        >>> errors["mean_error"] < 0.01  # Verify accuracy
+        True
+        >>> errors["max_error"] < 0.02
+        True
+        """
+        vertices = self.get_vertices()
+        if len(vertices) == 0:
+            return {"mean_error": 0.0, "max_error": 0.0}
+
+        # Sample vertices if there are many
+        if len(vertices) > n_samples:
+            indices = np.random.choice(len(vertices), size=n_samples, replace=False)
+            vertices = vertices[indices]
+
+        # Evaluate function at vertices
+        values = self.func(vertices[:, 0], vertices[:, 1], vertices[:, 2])
+        errors = np.abs(values - self.level)
+
+        return {
+            "mean_error": float(np.mean(errors)),
+            "max_error": float(np.max(errors)),
+        }
diff --git a/pyproject.toml b/pyproject.toml
index a0cfc3d11f..fac129ae0c 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -38,6 +38,7 @@ dependencies = [
     "numpy>=2.1 ; python_full_version >= '3.10'",
     "pillow>=9.1",
     "pycairo>=1.13,<2.0.0",
+    "pymcubes>=0.1.0",
     "pydub>=0.20.0",
     "pygments>=2.0.0",
     "rich>=12.0.0",
diff --git a/tests/module/mobject/test_implicit_surface.py b/tests/module/mobject/test_implicit_surface.py
new file mode 100644
index 0000000000..8a13477927
--- /dev/null
+++ b/tests/module/mobject/test_implicit_surface.py
@@ -0,0 +1,196 @@
+"""Unit tests for ImplicitSurface."""
+
+from __future__ import annotations
+
+import numpy as np
+
+from manim import GREEN, ImplicitSurface
+
+
+class TestImplicitSurface:
+    """Tests for the ImplicitSurface class."""
+
+    def test_sphere_creation(self):
+        """Test that a sphere can be created."""
+
+        def sphere_func(x, y, z):
+            return x**2 + y**2 + z**2 - 1.0
+
+        surface = ImplicitSurface(
+            sphere_func,
+            x_range=(-1.5, 1.5),
+            y_range=(-1.5, 1.5),
+            z_range=(-1.5, 1.5),
+            resolution=32,
+        )
+
+        # Should have triangular faces
+        assert len(surface.submobjects) > 0
+
+    def test_sphere_vertices_on_surface(self):
+        """Test that sphere vertices approximately satisfy the implicit equation."""
+
+        def sphere_func(x, y, z):
+            return x**2 + y**2 + z**2 - 1.0
+
+        surface = ImplicitSurface(
+            sphere_func,
+            x_range=(-1.5, 1.5),
+            y_range=(-1.5, 1.5),
+            z_range=(-1.5, 1.5),
+            resolution=40,
+        )
+
+        errors = surface.verify_surface()
+        # Mean error should be small
+        assert errors["mean_error"] < 0.1
+        # Max error should be reasonable
+        assert errors["max_error"] < 0.2
+
+    def test_no_surface_in_domain(self):
+        """Test that no triangles are created when surface is outside domain."""
+
+        def f(x, y, z):
+            return x**2 + y**2 + z**2 - 100.0  # Sphere of radius 10
+
+        surface = ImplicitSurface(
+            f,
+            x_range=(-1, 1),
+            y_range=(-1, 1),
+            z_range=(-1, 1),
+            resolution=20,
+        )
+
+        # No triangles should be created
+        assert len(surface.submobjects) == 0
+
+    def test_level_parameter(self):
+        """Test that the level parameter works correctly."""
+
+        def f(x, y, z):
+            return x**2 + y**2 + z**2
+
+        # Create sphere of radius 2 using level=4
+        surface = ImplicitSurface(
+            f,
+            x_range=(-2.5, 2.5),
+            y_range=(-2.5, 2.5),
+            z_range=(-2.5, 2.5),
+            level=4.0,  # radius squared = 4, so radius = 2
+            resolution=32,
+        )
+
+        assert len(surface.submobjects) > 0
+
+        # Check vertices have approximately radius 2
+        vertices = surface.get_vertices()
+        if len(vertices) > 0:
+            radii = np.sqrt(np.sum(vertices**2, axis=1))
+            mean_radius = np.mean(radii)
+            assert abs(mean_radius - 2.0) < 0.2
+
+    def test_color_parameters(self):
+        """Test that color parameters are applied."""
+
+        def f(x, y, z):
+            return x**2 + y**2 + z**2 - 1.0
+
+        surface = ImplicitSurface(
+            f,
+            resolution=20,
+            fill_color=GREEN,
+            fill_opacity=0.5,
+            stroke_width=1.0,
+        )
+
+        if len(surface.submobjects) > 0:
+            # Check that faces have the expected properties
+            face = surface.submobjects[0]
+            assert face.fill_opacity == 0.5
+            assert face.stroke_width == 1.0
+
+    def test_cylinder_surface(self):
+        """Test creating a cylinder using implicit function."""
+
+        def cylinder_func(x, y, z):
+            return x**2 + y**2 - 1.0
+
+        surface = ImplicitSurface(
+            cylinder_func,
+            x_range=(-1.5, 1.5),
+            y_range=(-1.5, 1.5),
+            z_range=(-1, 1),
+            resolution=32,
+        )
+
+        assert len(surface.submobjects) > 0
+
+    def test_torus_surface(self):
+        """Test creating a torus using implicit function."""
+        R, r = 2.0, 0.5
+
+        def torus_func(x, y, z):
+            return (np.sqrt(x**2 + y**2) - R) ** 2 + z**2 - r**2
+
+        surface = ImplicitSurface(
+            torus_func,
+            x_range=(-3, 3),
+            y_range=(-3, 3),
+            z_range=(-1, 1),
+            resolution=40,
+        )
+
+        assert len(surface.submobjects) > 0
+
+    def test_get_vertices(self):
+        """Test getting vertices from the surface."""
+
+        def f(x, y, z):
+            return x**2 + y**2 + z**2 - 1.0
+
+        surface = ImplicitSurface(f, resolution=32)
+        vertices = surface.get_vertices()
+
+        # Should have vertices
+        assert len(vertices) > 0
+        # Each vertex should be 3D
+        assert vertices.shape[1] == 3
+
+    def test_empty_surface_get_vertices(self):
+        """Test getting vertices from an empty surface."""
+
+        def f(x, y, z):
+            return x**2 + y**2 + z**2 - 100.0
+
+        surface = ImplicitSurface(
+            f,
+            x_range=(-1, 1),
+            y_range=(-1, 1),
+            z_range=(-1, 1),
+            resolution=20,
+        )
+
+        vertices = surface.get_vertices()
+        assert len(vertices) == 0
+
+    def test_plane_surface(self):
+        """Test creating a plane using implicit function."""
+
+        def plane_func(x, y, z):
+            return z - 0.5  # Plane at z = 0.5
+
+        surface = ImplicitSurface(
+            plane_func,
+            x_range=(-1, 1),
+            y_range=(-1, 1),
+            z_range=(0, 1),
+            resolution=32,
+        )
+
+        assert len(surface.submobjects) > 0
+
+        # Vertices should have z close to 0.5
+        vertices = surface.get_vertices()
+        if len(vertices) > 0:
+            z_values = vertices[:, 2]
+            assert np.allclose(z_values, 0.5, atol=0.1)
diff --git a/uv.lock b/uv.lock
index caad3683c4..89fa298b02 100644
--- a/uv.lock
+++ b/uv.lock
@@ -1,5 +1,5 @@
 version = 1
-revision = 3
+revision = 2
 requires-python = ">=3.10"
 resolution-markers = [
     "python_full_version >= '3.13'",
@@ -1350,6 +1350,7 @@ dependencies = [
     { name = "pycairo" },
     { name = "pydub" },
     { name = "pygments" },
+    { name = "pymcubes" },
     { name = "rich" },
     { name = "scipy", version = "1.13.1", source = { registry = "https://pypi.org/simple" }, marker = "python_full_version < '3.13'" },
     { name = "scipy", version = "1.15.1", source = { registry = "https://pypi.org/simple" }, marker = "python_full_version >= '3.13'" },
@@ -1419,6 +1420,7 @@ requires-dist = [
     { name = "pycairo", specifier = ">=1.13,<2.0.0" },
     { name = "pydub", specifier = ">=0.20.0" },
     { name = "pygments", specifier = ">=2.0.0" },
+    { name = "pymcubes", specifier = ">=0.1.0" },
     { name = "rich", specifier = ">=12.0.0" },
     { name = "scipy", specifier = ">=1.13.0" },
     { name = "scipy", marker = "python_full_version >= '3.13'", specifier = ">=1.14.0" },
@@ -2299,6 +2301,31 @@ crypto = [
     { name = "cryptography" },
 ]
 
+[[package]]
+name = "pymcubes"
+version = "0.1.6"
+source = { registry = "https://pypi.org/simple" }
+dependencies = [
+    { name = "numpy" },
+    { name = "scipy", version = "1.13.1", source = { registry = "https://pypi.org/simple" }, marker = "python_full_version < '3.13'" },
+    { name = "scipy", version = "1.15.1", source = { registry = "https://pypi.org/simple" }, marker = "python_full_version >= '3.13'" },
+]
+sdist = { url = "https://files.pythonhosted.org/packages/df/fb/8660bc09b05e5b7209377a476e937f36c2e4353a5300e6c71283fd6f8c2b/pymcubes-0.1.6.tar.gz", hash = "sha256:1a72476a22a4c98907acec0c5ef8d46f4451fa987c7efeaf970990775826ed8b", size = 109395, upload-time = "2024-07-03T23:14:43.434Z" }
+wheels = [
+    { url = "https://files.pythonhosted.org/packages/e9/40/cd07da8d67a3e70a22e3fcd86280a287dc95f9118b9935c7fb6771850be2/PyMCubes-0.1.6-cp310-cp310-macosx_10_9_x86_64.whl", hash = "sha256:8b66c714abb554dbb0cc3fa75fc1fd24e27dc2cc618a2685832b825156848792", size = 54392, upload-time = "2024-07-03T23:17:13.756Z" },
+    { url = "https://files.pythonhosted.org/packages/4e/ef/0ba7253b8ea48bc6207ef7f5dbcf61df7237bf8ae0c1bb86ce39c6432cfa/PyMCubes-0.1.6-cp310-cp310-macosx_11_0_arm64.whl", hash = "sha256:e37a81ffff099e8d275c4c973f34978934720120a02d06198c410f48c9f6b7d0", size = 52470, upload-time = "2024-07-03T23:15:24.635Z" },
+    { url = "https://files.pythonhosted.org/packages/5e/46/c4d6d8faed01bb63be83570a9d8d4b88963d82d84085fd95cf5423319d23/PyMCubes-0.1.6-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:16ed358f20f514def6b4669d4714f86429c609421f167a2fd28731b71adf381a", size = 318360, upload-time = "2024-07-03T23:15:52.994Z" },
+    { url = "https://files.pythonhosted.org/packages/60/f9/db2002e639681476d863674e1f52290f02d7a701700676f9f53c5f8ea4d2/PyMCubes-0.1.6-cp310-cp310-win_amd64.whl", hash = "sha256:8433ecf1ef8dca963056d0b5e8ce9d6a1d5160764a199fd1a375a4ab58e89a3c", size = 50022, upload-time = "2024-07-03T23:18:51.281Z" },
+    { url = "https://files.pythonhosted.org/packages/62/d3/4eb02520fd20558ec581af2bd5f670fd58c31c5630d90a7e30f1be501222/PyMCubes-0.1.6-cp311-cp311-macosx_10_9_x86_64.whl", hash = "sha256:e859a7aaa97bfc909079ce1ecbc12f0bc61b49b8b01ee0d59eaabf9a7f266a9f", size = 54338, upload-time = "2024-07-03T23:17:15.464Z" },
+    { url = "https://files.pythonhosted.org/packages/3a/e8/c39c7eedb7a2b514092ea8017b11f1c3ad3cc9398e5afa7568739a859447/PyMCubes-0.1.6-cp311-cp311-macosx_11_0_arm64.whl", hash = "sha256:8f1246d4ef9dee648a111060aac95f0b99f9b36c937436e75af4cae57b5f50f7", size = 52324, upload-time = "2024-07-03T23:15:26.161Z" },
+    { url = "https://files.pythonhosted.org/packages/2a/e6/5796cf17a909d6048c0c2cfd70577029fb39c15e5aa098df5336b316a6c7/PyMCubes-0.1.6-cp311-cp311-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:323b1d6454e3092ce2ba42c0ce69a753cb646fbda33d5274f03514a68f4441c8", size = 336818, upload-time = "2024-07-03T23:15:54.554Z" },
+    { url = "https://files.pythonhosted.org/packages/bd/4c/30cfeb335589ddcb72db0f4f68e5fce703ada8d05e4c2b37b5125bbae1a9/PyMCubes-0.1.6-cp311-cp311-win_amd64.whl", hash = "sha256:4e326cba1471047567e016c39fa2b4934dd451e9f1feec63a35e23bbcab1d34c", size = 50006, upload-time = "2024-07-03T23:18:52.403Z" },
+    { url = "https://files.pythonhosted.org/packages/27/17/71df57cadc0f34d31f0348f6067dea19de3c100e9aed069b792e8c4ca8fe/PyMCubes-0.1.6-cp312-cp312-macosx_10_9_x86_64.whl", hash = "sha256:6dea0bfb68778d7f9becd87c86a396c9ffbb8f322c184967fa2a48d430974513", size = 53794, upload-time = "2024-07-03T23:17:16.767Z" },
+    { url = "https://files.pythonhosted.org/packages/7f/50/78eeda22066c42c6a5c77b7b762e670880b8ec970b1f38ba11d575c3ba31/PyMCubes-0.1.6-cp312-cp312-macosx_11_0_arm64.whl", hash = "sha256:138eaa87d2b8daf7488102070191e8e99254d48b6f08577a618db6e95fd53f90", size = 51547, upload-time = "2024-07-03T23:15:28.124Z" },
+    { url = "https://files.pythonhosted.org/packages/b5/13/63d018ea286c95b758f8b5cbe97d88a16b38f5672411138287bd2123d877/PyMCubes-0.1.6-cp312-cp312-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:ea366a2064af0846093e0ad3f9035e375f4b14b62bb565c95dcc8dcaf78308a5", size = 339385, upload-time = "2024-07-03T23:15:56.087Z" },
+    { url = "https://files.pythonhosted.org/packages/f6/3e/a9a9f624a364166b1b4513f5ae8627fb59eea64eccf3e5a91c82c2b7b44e/PyMCubes-0.1.6-cp312-cp312-win_amd64.whl", hash = "sha256:8926475690ce1e0be2460c941a51cdda6eed6f10986073d23fc67535bc70aa78", size = 49968, upload-time = "2024-07-03T23:18:53.376Z" },
+]
+
 [[package]]
 name = "pynacl"
 version = "1.5.0"