PR: Add meshgrid

spec/API_specification/creation_functions.md

@@ -289,6 +289,49 @@ Returns evenly spaced numbers over a specified interval.
 
     -   a one-dimensional array containing evenly spaced values.
 
+(function-meshgrid)=
+### meshgrid(*arrays, /, indexing='xy')
+
+Returns coordinate matrices from coordinate vectors.
+
+### Special cases
+
+-    For `n` one dimensional arrays `x1, x2, ..., xn` with lengths `Ni = len(xi)`, this function returns `(N1, N2, N3, ..., Nn)` shaped arrays if the indexing is 'ij' or `(N2, N1, N3, ..., Nn)` shaped arrays if indexing 'xy'.
+
+-    The 0-D and 1-D case, the indexing keyword has no effect.
+
+### Parameters
+
+-    **\*arrays**: _<array>_
+
+     -   One dimensional arrays representing the coordinates of the grid.
+
+-    **indexing**: _Optional\[ str ]_
+
+     -   Cartesian 'xy' by default or 'ij' for matrix indexing of the output.
+
+### Returns
+
+-    **out**: _Tuple\[ <array>, ... ]_
+
+      -    List of N arrays with rank `max(2, N)`.
+
+### Notes
+
+This function supports indexing conventions using the `indexing` keyword. In the case of `xy` the function will return a meshgrid with cartesian indexing, while `ij` will use matrix indexing. The difference is illustrated by the following code snippet:
+
+```
+xv, yv = meshgrid(x, y, indexing='xy')
+for i in range(nx):
+    for j in range(ny):
+        # treat xv[j, i], yv[j, i]
+
+xv, yv = meshgrid(x, y, indexing='ij')
+for i in range(nx):
+    for j in range(ny):
+        # treat xv[i, j], yv[i, j]
+```
+
 (function-ones)=
 ### ones(shape, /, *, dtype=None, device=None)