diff --git a/doc/Programs/NumpyCodes/meshgrid.py b/doc/Programs/NumpyCodes/meshgrid.py
new file mode 100644
index 000000000..c21d4c08f
--- /dev/null
+++ b/doc/Programs/NumpyCodes/meshgrid.py
@@ -0,0 +1,38 @@
+import numpy as np
+from random import random, seed
+
+
+# Make data.
+x = np.arange(0, 1, 0.5)
+y = np.arange(0, 1, 0.5)
+x, y = np.meshgrid(x,y)
+
+def FrankeFunction(x,y):
+    term1 = 0.75*np.exp(-(0.25*(9*x-2)**2) - 0.25*((9*y-2)**2))
+    term2 = 0.75*np.exp(-((9*x+1)**2)/49.0 - 0.1*(9*y+1))
+    term3 = 0.5*np.exp(-(9*x-7)**2/4.0 - 0.25*((9*y-3)**2))
+    term4 = -0.2*np.exp(-(9*x-4)**2 - (9*y-7)**2)
+    return term1 + term2 + term3 + term4
+
+
+def create_X(x, y, n ):
+	if len(x.shape) > 1:
+		x = np.ravel(x)
+		y = np.ravel(y)
+
+	N = len(x)
+	l = int((n+1)*(n+2)/2)		# Number of elements in beta
+	X = np.ones((N,l))
+
+	for i in range(1,n+1):
+		q = int((i)*(i+1)/2)
+		for k in range(i+1):
+			X[:,q+k] = (x**(i-k))*(y**k)
+
+	return X
+
+# Making meshgrid of datapoints and compute Franke's function
+p = 2
+z = FrankeFunction(x, y)
+X = create_X(x, y, p)        
+print(X)