diff --git a/C/main.c b/C/main.c
index e1a7fcfc15cfc8577fb5902837783f2c0fcc55ee..478f84ccfef97f3d3e827cac348d3c176465897f 100644
--- a/C/main.c
+++ b/C/main.c
@@ -10,23 +10,70 @@
  *  @param x A double variable.
  *  @return f(x) = 2.0*sin(x) - 3.0*cos(x)
  */
-double my_function(double x)
+double my_1d_function(double x)
 {
     return 2.0 * sin(x) - 3.0 * cos(x);
 }
 
-int main()
+/** @brief 
+ *  An exemple of mathematical function 
+ *  in 2 dimension.
+ * 
+ *  @param x The first argument.
+ *  @param y The second argument.
+ *  @return f(x,y) = x^2 + y^2
+ */
+double my_2d_function(double x, double y)
+{
+    return x * x + y * y;
+}
+
+/** @brief 
+ *  Generate an example of 1D data,
+ *  wich is stored into two vec files.
+ */
+void generate_1D_data()
 {
     // Create a vector X = [0,1,2...99]
     double_vector_t *X = iota(100);
-    // Create a vector Y = my_function(x)
-    double_vector_t *Y = apply_function(X, my_function);
+    // Create a vector Y = my_1d_function(x)
+    double_vector_t *Y = apply_1d_function(X, my_1d_function);
 
     // Export our vectors into files
-    export_vector("../X.vec", X);
-    export_vector("../Y.vec", Y);
+    export_vector("../X_1D.vec", X);
+    export_vector("../Y_1D.vec", Y);
 
     // Free our vectors
     destroy_vector(&Y);
     destroy_vector(&X);
+}
+
+/** @brief 
+ *  Generate an example of 2D data,
+ *  wich is stored into three vec files.
+ */
+void generate_2D_data()
+{
+    // Create a vector X = [0,1,2...99]
+    double_vector_t *X = iota(100);
+    // Create a vector Y = [0,1,2...99]
+    double_vector_t *Y = iota(100);
+    // Create a vector Z = my_2d_function(X, Y)
+    double_vector_t *Z = apply_2d_function(X, Y, my_2d_function);
+
+    // Export our vectors into files
+    export_vector("../X_2D.vec", X);
+    export_vector("../Y_2D.vec", Y);
+    export_vector("../Z_2D.vec", Z);
+
+    // Free our vectors
+    destroy_vector(&Z);
+    destroy_vector(&Y);
+    destroy_vector(&X);
+}
+
+int main()
+{
+    generate_1D_data();
+    generate_2D_data();
 }
\ No newline at end of file
diff --git a/C/vector.c b/C/vector.c
index 431f1cf4b067ab9cce41239db09bf18501c0af56..31afaddba5bd2eb45f37c981e815965b725ca08c 100644
--- a/C/vector.c
+++ b/C/vector.c
@@ -56,18 +56,43 @@ double_vector_t *iota(uint32_t N)
  *  to a given vector, and return the
  *  result in a new vector.
  * 
- *  @param vec The argument vector
- *  @param f   The 1d function to apply
+ *  @param X The argument vector
+ *  @param f The 1d function to apply
  *  @return A dynamically allocated vector : f(X)
  */
-double_vector_t *apply_function(double_vector_t *vec, double_function_t f)
+double_vector_t *apply_1d_function(double_vector_t *X, double_1d_function_t f)
 {
-    double_vector_t *res = init_vector(vec->N);
-    for (uint32_t i = 0; i < vec->N; i++)
+    double_vector_t *Y = init_vector(X->N);
+    for (uint32_t i = 0; i < X->N; i++)
     {
-        res->components[i] = f(vec->components[i]);
+        Y->components[i] = f(X->components[i]);
     }
-    return res;
+    return Y;
+}
+/** @brief 
+ *  Apply a 2d function element-wise
+ *  to given vectors, and return the
+ *  result in a new vector.
+ * 
+ *  @param X The first argument vector
+ *  @param Y The second argument vector
+ *  @param f The 2d function to apply
+ *  @return A dynamically allocated vector : f(X,Y)
+ */
+double_vector_t *apply_2d_function(double_vector_t *X, double_vector_t *Y, double_2d_function_t f)
+{
+    if (X->N != Y->N)
+    {
+        fprintf(stderr, "Bad shape len(X) = %d and len(Y) = %d\n", X->N, Y->N);
+        return NULL;
+    }
+
+    double_vector_t *Z = init_vector(X->N);
+    for (uint32_t i = 0; i < X->N; i++)
+    {
+        Z->components[i] = f(X->components[i], Y->components[i]);
+    }
+    return Z;
 }
 /** @brief 
  *  Export a vector into a file.
@@ -78,6 +103,11 @@ double_vector_t *apply_function(double_vector_t *vec, double_function_t f)
 void export_vector(const char *filename, double_vector_t *vec)
 {
     FILE *output = fopen(filename, "w");
+    if (output == NULL)
+    {
+        fprintf(stderr, "Can't open output file");
+        exit(EXIT_FAILURE);
+    }
 
     vector_metadata_t metadata;
     metadata.endianness = get_endianness();
diff --git a/C/vector.h b/C/vector.h
index 4b391f4581b2d681778177bc6b056d64674b6c16..d1df6114d5aef1f4d67be977e003228f76f312c1 100644
--- a/C/vector.h
+++ b/C/vector.h
@@ -6,7 +6,9 @@ typedef struct double_vector
     double *components;
 } double_vector_t;
 // Function pointer, example : double f(double x);
-typedef double (*double_function_t)(double);
+typedef double (*double_1d_function_t)(double);
+// Function pointer, example : double f(double x, double y);
+typedef double (*double_2d_function_t)(double, double);
 
 /*
 * The attribute "packed" tells the compiler,
@@ -42,11 +44,22 @@ double_vector_t *iota(uint32_t N);
  *  to a given vector, and return the
  *  result in a new vector.
  * 
- *  @param vec The argument vector
- *  @param f   The 1d function to apply
+ *  @param X The argument vector
+ *  @param f The 1d function to apply
  *  @return A dynamically allocated vector : f(X)
  */
-double_vector_t *apply_function(double_vector_t *vec, double_function_t f);
+double_vector_t *apply_1d_function(double_vector_t *X, double_1d_function_t f);
+/** @brief 
+ *  Apply a 2d function element-wise
+ *  to given vectors, and return the
+ *  result in a new vector.
+ * 
+ *  @param X The first argument vector
+ *  @param Y The second argument vector
+ *  @param f The 2d function to apply
+ *  @return A dynamically allocated vector : f(X,Y)
+ */
+double_vector_t *apply_2d_function(double_vector_t *X, double_vector_t *Y, double_2d_function_t f);
 /** @brief 
  *  Export a vector into a file.
  * 
diff --git a/Python/example.py b/Python/example.py
index e40a55698493ef70de538cc86c46d2f79af95e43..f80b7ba77892f344639ca4cb9662f429e10eeddc 100644
--- a/Python/example.py
+++ b/Python/example.py
@@ -1,19 +1,51 @@
+from mpl_toolkits.mplot3d import Axes3D
+from matplotlib import cm
 from matplotlib import pyplot as plt
 from load_vec import load_vector
+import numpy as np
 
-X = load_vector("../X.vec")
-Y = load_vector("../Y.vec")
 
-type_of_data = 'curve'
+TYPE_OF_DATA = 'curve'
 
-if type_of_data == 'curve':
-    plt.plot(X, Y, label="my curve")
-else:
-    plt.scatter(X, Y, marker='x', label="my points")
 
-plt.title("My data")
-plt.xlabel("X")
-plt.ylabel("Y")
-plt.legend(loc="upper right")
+def show_1d():
+    X = load_vector("../X_1D.vec")
+    Y = load_vector("../Y_1D.vec")
 
-plt.show()
+    if TYPE_OF_DATA == 'curve':
+        plt.plot(X, Y, c='g', label="my curve")
+    else:
+        plt.scatter(X, Y, c='r', marker='x', label="my points")
+
+    plt.title("My data")
+    plt.xlabel("X")
+    plt.ylabel("Y")
+    plt.legend(loc="upper right")
+
+    plt.show()
+
+
+def show_2d():
+    fig = plt.figure()
+    ax = fig.gca(projection='3d')
+
+    X = load_vector("../X_2D.vec")
+    Y = load_vector("../Y_2D.vec")
+    Z = load_vector("../Z_2D.vec")
+
+    if TYPE_OF_DATA == 'curve':
+        ax.plot(X, Y, Z, c='r', label="my curve")
+    else:
+        ax.scatter(X, Y, Z, c='r', marker='x', label="my points")
+
+    ax.set_title("My data")
+    ax.set_xlabel("X")
+    ax.set_ylabel("Y")
+    ax.set_zlabel("Z")
+    ax.legend(loc="upper right")
+
+    plt.show()
+
+
+show_1d()
+show_2d()