phpmlkit · CodeWithKyrian · May 30, 2026 · May 30, 2026
diff --git a/SPEC.md b/SPEC.md
@@ -73,6 +73,7 @@ The primary array type representing N-dimensional arrays.
 - [x] 3.2.2: `NDArray::linspace($start, $stop, $num = 50, $endpoint = true, $dtype = null)` - Linear spacing
 - [x] 3.2.3: `NDArray::logspace($start, $stop, $num = 50, $base = 10.0, $dtype = null)` - Logarithmic spacing
 - [x] 3.2.4: `NDArray::geomspace($start, $stop, $num = 50, $dtype = null)` - Geometric spacing
+- [x] 3.2.5: `NDArray::meshgrid($arrays, $indexing = 'xy', $sparse = false)` - Coordinate matrices from coordinate vectors
 
 ### 3.3 Random Array Creation (REQ-3.3)
 **Priority**: HIGH

diff --git a/docs/api/array-creation.md b/docs/api/array-creation.md
@@ -453,6 +453,81 @@ echo $floats;  // [0. 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9]
 
 ---
 
+## NDArray::meshgrid()
+
+Create coordinate matrices from one-dimensional coordinate vectors.
+
+```php
+public static function meshgrid(array $arrays, string $indexing = 'xy', bool $sparse = false): array
+```
+
+Given one or more coordinate vectors, `meshgrid()` returns one grid for each input. Dense grids are expanded to the full coordinate shape. Sparse grids keep singleton dimensions and avoid repeating coordinate values.
+
+**Parameters:**
+- `array $arrays` - One-dimensional coordinate vectors. Each entry may be an `NDArray` or a PHP array.
+- `string $indexing` - Coordinate indexing mode: `'xy'` for Cartesian indexing or `'ij'` for matrix indexing. Default: `'xy'`.
+- `bool $sparse` - Whether to return sparse grids. Default: `false`.
+
+**Returns:** Array of NDArrays, one coordinate grid for each input vector. Each output keeps the dtype of its corresponding input vector.
+
+**Throws:**
+- `ShapeException` - If no input arrays are provided.
+- `ShapeException` - If `indexing` is not `'xy'` or `'ij'`.
+- `ShapeException` - If any input is not one-dimensional.
+
+**Examples:**
+
+```php
+$x = NDArray::array([1, 2, 3]);
+$y = NDArray::array([10, 20]);
+
+[$xx, $yy] = NDArray::meshgrid([$x, $y]);
+
+print_r($xx->shape());
+// Output: [2, 3]
+
+print_r($xx->toArray());
+// Output: [[1, 2, 3], [1, 2, 3]]
+
+print_r($yy->toArray());
+// Output: [[10, 10, 10], [20, 20, 20]]
+```
+
+Use matrix indexing when the first output axis should correspond to the first input vector:
+
+```php
+[$xx, $yy] = NDArray::meshgrid([$x, $y], indexing: 'ij');
+
+print_r($xx->shape());
+// Output: [3, 2]
+
+print_r($xx->toArray());
+// Output: [[1, 1], [2, 2], [3, 3]]
+
+print_r($yy->toArray());
+// Output: [[10, 20], [10, 20], [10, 20]]
+```
+
+Sparse grids keep only the dimensions needed for each coordinate vector:
+
+```php
+[$xx, $yy] = NDArray::meshgrid([$x, $y], sparse: true);
+
+print_r($xx->shape());
+// Output: [1, 3]
+
+print_r($yy->shape());
+// Output: [2, 1]
+```
+
+**See Also:**
+- [arange()](#ndarray-arange)
+- [linspace()](#ndarray-linspace)
+- [reshape()](/api/array-manipulation#reshape)
+- [tile()](/api/array-manipulation#tile)
+
+---
+
 ## NDArray::linspace()
 
 Create linearly spaced values.

diff --git a/docs/api/global-functions.md b/docs/api/global-functions.md
@@ -56,6 +56,7 @@ You can combine imports from the base namespace and sub-namespaces in one `use f
 | `linspace` | `NDArray::linspace()` | [Array Creation — linspace](/api/array-creation#ndarray-linspace) |
 | `logspace` | `NDArray::logspace()` | [Array Creation - logspace](/api/array-creation#ndarray-logspace) |
 | `geomspace` | `NDArray::geomspace()` | [Array Creation - geomspace](/api/array-creation#ndarray-geomspace) |
+| `meshgrid` | `NDArray::meshgrid()` | [Array Creation - meshgrid](/api/array-creation#ndarray-meshgrid) |
 | `random` | `NDArray::random()` | [Array Creation - random](/api/array-creation#ndarray-random) |
 | `random_int` | `NDArray::randomInt()` | [Array Creation](/api/array-creation) |
 | `randn` | `NDArray::randn()` | [Array Creation](/api/array-creation) |

diff --git a/docs/guide/fundamentals/operations.md b/docs/guide/fundamentals/operations.md
@@ -285,6 +285,11 @@ $matrix->ravel();    // View if possible
 $arr = NDArray::array([1, 2, 3]);
 $arr->expandDims(0);   // Shape [1, 3]
 $arr->squeeze();        // Remove size-1 dimensions
+
+// Coordinate grids
+$x = NDArray::array([1, 2, 3]);
+$y = NDArray::array([10, 20]);
+[$xx, $yy] = NDArray::meshgrid([$x, $y]);  // Shape [2, 3] for both grids
 ```
 
 **Views vs Copies:**

diff --git a/docs/guide/getting-started/quick-start.md b/docs/guide/getting-started/quick-start.md
@@ -79,6 +79,11 @@ $step = NDArray::arange(0, 10, 2); // [0 2 4 6 8]
 // linspace: linear spacing
 $linear = NDArray::linspace(0, 1, 5);  // [0. 0.25 0.5 0.75 1.]
 
+// Coordinate grids
+$x = NDArray::array([1, 2, 3]);
+$y = NDArray::array([10, 20]);
+[$xx, $yy] = NDArray::meshgrid([$x, $y]);  // Both grids have shape [2, 3]
+
 // Random arrays
 $random = NDArray::random([3, 3]);         // Uniform [0, 1)
 $normal = NDArray::randn([3, 3]);          // Standard normal distribution

diff --git a/src/Functions.php b/src/Functions.php
@@ -128,6 +128,20 @@ function full_like(NDArray $array, bool|Complex|float|int $value, ?DType $dtype
         return NDArray::fullLike($array, $value, $dtype);
     }
 
+    /**
+     * Create coordinate matrices from one-dimensional coordinate vectors.
+     *
+     * @param array<array<mixed>|NDArray> $arrays   One-dimensional coordinate vectors
+     * @param string                      $indexing Cartesian ('xy') or matrix ('ij') indexing
+     * @param bool                        $sparse   Whether to return sparse coordinate grids
+     *
+     * @return array<NDArray> Coordinate grids, one for each input vector
+     */
+    function meshgrid(array $arrays, string $indexing = 'xy', bool $sparse = false): array
+    {
+        return NDArray::meshgrid($arrays, $indexing, $sparse);
+    }
+
     /**
      * Create a 2D identity matrix.
      *

diff --git a/src/Traits/CreatesArrays.php b/src/Traits/CreatesArrays.php
@@ -296,6 +296,92 @@ public static function fullLike(self $array, bool|Complex|float|int $value, ?DTy
         return self::full($value, $array->shape(), $dtype);
     }
 
+    /**
+     * Create coordinate matrices from one-dimensional coordinate vectors.
+     *
+     * Dense grids are materialized with tile(). Sparse grids keep singleton dimensions
+     * and avoid expanding repeated coordinate values.
+     *
+     * @param array<array<mixed>|self> $arrays   One-dimensional coordinate vectors
+     * @param 'ij'|'xy'                $indexing Cartesian ('xy') or matrix ('ij') indexing
+     * @param bool                     $sparse   Whether to return sparse coordinate grids
+     *
+     * @return array<self> Coordinate grids, one for each input vector
+     *
+     * @throws ShapeException If no arrays are provided, an input is not one-dimensional, or indexing is invalid
+     */
+    public static function meshgrid(array $arrays, string $indexing = 'xy', bool $sparse = false): array
+    {
+        if ([] === $arrays) {
+            throw new ShapeException('meshgrid() requires at least one input array');
+        }
+
+        if (!\in_array($indexing, ['xy', 'ij'], true)) {
+            throw new ShapeException("meshgrid() indexing must be 'xy' or 'ij', got '{$indexing}'");
+        }
+
+        $vectors = [];
+        foreach ($arrays as $i => $array) {
+            if ($array instanceof self) {
+                $vector = $array;
+            } elseif (\is_array($array)) {
+                $vector = self::array($array);
+            } else {
+                throw new ShapeException(
+                    'meshgrid() inputs must be NDArrays or PHP arrays, got '.get_debug_type($array).' at index '.$i
+                );
+            }
+
+            if (1 !== $vector->ndim()) {
+                throw new ShapeException(
+                    "meshgrid() inputs must be one-dimensional, input {$i} has {$vector->ndim()} dimensions"
+                );
+            }
+
+            $vectors[] = $vector;
+        }
+
+        $count = \count($vectors);
+        $lengths = array_map(static fn (self $vector): int => $vector->size(), $vectors);
+        $outputShape = $lengths;
+
+        if ('xy' === $indexing && $count >= 2) {
+            $outputShape = $lengths;
+            $outputShape[0] = $lengths[1];
+            $outputShape[1] = $lengths[0];
+        }
+
+        $grids = [];
+        foreach ($vectors as $i => $vector) {
+            $axis = $i;
+            if ('xy' === $indexing && $count >= 2) {
+                $axis = match ($i) {
+                    0 => 1,
+                    1 => 0,
+                    default => $i,
+                };
+            }
+
+            $baseShape = array_fill(0, $count, 1);
+            $baseShape[$axis] = $lengths[$i];
+
+            $grid = $vector->reshape($baseShape);
+
+            if (!$sparse) {
+                $reps = [];
+                foreach ($outputShape as $dim => $dimSize) {
+                    $reps[] = $dim === $axis ? 1 : $dimSize;
+                }
+
+                $grid = $grid->tile($reps);
+            }
+
+            $grids[] = $grid;
+        }
+
+        return $grids;
+    }
+
     /**
      * Create a 2D identity matrix.
      *

diff --git a/tests/Unit/CreationTest.php b/tests/Unit/CreationTest.php
@@ -10,6 +10,8 @@
 use PhpMlKit\NDArray\NDArray;
 use PHPUnit\Framework\TestCase;
 
+use function PhpMlKit\NDArray\meshgrid;
+
 /**
  * @internal
  *
@@ -69,6 +71,114 @@ public function testFullInferred(): void
         $this->assertSame([true, true], $arrBool->toArray());
     }
 
+    public function testMeshgridDenseXyIndexing(): void
+    {
+        $x = NDArray::array([1, 2, 3], DType::Int64);
+        $y = NDArray::array([10, 20], DType::Int64);
+
+        [$xx, $yy] = NDArray::meshgrid([$x, $y]);
+
+        $this->assertSame([2, 3], $xx->shape());
+        $this->assertSame([2, 3], $yy->shape());
+        $this->assertSame([[1, 2, 3], [1, 2, 3]], $xx->toArray());
+        $this->assertSame([[10, 10, 10], [20, 20, 20]], $yy->toArray());
+    }
+
+    public function testMeshgridDenseIjIndexing(): void
+    {
+        $x = NDArray::array([1, 2, 3], DType::Int64);
+        $y = NDArray::array([10, 20], DType::Int64);
+
+        [$xx, $yy] = NDArray::meshgrid([$x, $y], indexing: 'ij');
+
+        $this->assertSame([3, 2], $xx->shape());
+        $this->assertSame([3, 2], $yy->shape());
+        $this->assertSame([[1, 1], [2, 2], [3, 3]], $xx->toArray());
+        $this->assertSame([[10, 20], [10, 20], [10, 20]], $yy->toArray());
+    }
+
+    public function testMeshgridSparseXyIndexing(): void
+    {
+        $x = NDArray::array([1, 2, 3], DType::Int64);
+        $y = NDArray::array([10, 20], DType::Int64);
+
+        [$xx, $yy] = NDArray::meshgrid([$x, $y], sparse: true);
+
+        $this->assertSame([1, 3], $xx->shape());
+        $this->assertSame([2, 1], $yy->shape());
+        $this->assertSame([[1, 2, 3]], $xx->toArray());
+        $this->assertSame([[10], [20]], $yy->toArray());
+        $this->assertTrue($xx->isView());
+        $this->assertTrue($yy->isView());
+    }
+
+    public function testMeshgridThreeInputsXyIndexing(): void
+    {
+        [$xx, $yy, $zz] = NDArray::meshgrid([
+            NDArray::array([1, 2], DType::Int64),
+            NDArray::array([10, 20, 30], DType::Int64),
+            NDArray::array([100, 200], DType::Int64),
+        ]);
+
+        $this->assertSame([3, 2, 2], $xx->shape());
+        $this->assertSame([3, 2, 2], $yy->shape());
+        $this->assertSame([3, 2, 2], $zz->shape());
+        $this->assertSame(1, $xx->get(0, 0, 0));
+        $this->assertSame(2, $xx->get(0, 1, 0));
+        $this->assertSame(20, $yy->get(1, 0, 0));
+        $this->assertSame(200, $zz->get(0, 0, 1));
+    }
+
+    public function testMeshgridAcceptsPhpArraysAndPreservesDtypes(): void
+    {
+        $x = NDArray::array([1, 2], DType::Int32);
+        $y = NDArray::array([1.5, 2.5], DType::Float64);
+
+        [$xx, $yy] = NDArray::meshgrid([$x, $y]);
+        [$fromPhpX, $fromPhpY] = NDArray::meshgrid([[1, 2], [3, 4]], indexing: 'ij');
+
+        $this->assertSame(DType::Int32, $xx->dtype());
+        $this->assertSame(DType::Float64, $yy->dtype());
+        $this->assertSame([[1, 1], [2, 2]], $fromPhpX->toArray());
+        $this->assertSame([[3, 4], [3, 4]], $fromPhpY->toArray());
+    }
+
+    public function testMeshgridGlobalFunction(): void
+    {
+        [$xx, $yy] = meshgrid([
+            NDArray::array([1, 2], DType::Int64),
+            NDArray::array([3, 4], DType::Int64),
+        ], indexing: 'ij');
+
+        $this->assertSame([[1, 1], [2, 2]], $xx->toArray());
+        $this->assertSame([[3, 4], [3, 4]], $yy->toArray());
+    }
+
+    public function testMeshgridRejectsEmptyInputs(): void
+    {
+        $this->expectException(ShapeException::class);
+        $this->expectExceptionMessage('at least one input array');
+
+        NDArray::meshgrid([]);
+    }
+
+    public function testMeshgridRejectsInvalidIndexing(): void
+    {
+        $this->expectException(ShapeException::class);
+        $this->expectExceptionMessage("indexing must be 'xy' or 'ij'");
+
+        // @phpstan-ignore argument.type (Intentionally passing an invalid indexing value)
+        NDArray::meshgrid([NDArray::array([1, 2])], indexing: 'bad');
+    }
+
+    public function testMeshgridRejectsNonVectorInputs(): void
+    {
+        $this->expectException(ShapeException::class);
+        $this->expectExceptionMessage('one-dimensional');
+
+        NDArray::meshgrid([NDArray::array([[1, 2], [3, 4]])]);
+    }
+
     public function testEye(): void
     {
         $arr = NDArray::eye(3);