diff --git a/examples/python/README.rst b/examples/python/README.rst
new file mode 100644
index 000000000..cbe2ff4a3
--- /dev/null
+++ b/examples/python/README.rst
@@ -0,0 +1,9 @@
+In order to use the Akantu examples using the python mode, one has to source
+the file 'akantu_test_environement.sh':
+
+> source AKANTU_BUILD_DIR/akantu_environement.sh
+
+or for an installed akantu version:
+
+> source AKANTU_INSTALL_DIR/share/akantuVERSION/akantu_environement.sh
+
diff --git a/examples/python/custom-material-dynamics/bar.geo b/examples/python/custom-material-dynamics/bar.geo
new file mode 100644
index 000000000..32f173f35
--- /dev/null
+++ b/examples/python/custom-material-dynamics/bar.geo
@@ -0,0 +1,45 @@
+// Mesh size
+h  = 0.05;
+
+h1 = h;
+h2 = h;
+
+// Dimensions of the bar
+Lx = 10;
+Ly = 1;
+
+// ------------------------------------------
+// Geometry
+// ------------------------------------------
+
+Point(101) = { 0.0, -Ly/2, 0.0, h1};
+Point(102) = { Lx,  -Ly/2, 0.0, h2};
+
+Point(103) = { Lx,  0., 0.0,  h2};
+Point(104) = { Lx,  Ly/2., 0.0,  h2};
+
+Point(105) = { 0.0, Ly/2., 0.0,  h1};
+Point(106) = { 0.0, 0., 0.0,  h1};
+
+Line(101) = {101,102};
+Line(102) = {102,103};
+Line(103) = {103,104};
+Line(104) = {104,105};
+Line(105) = {105,106};
+Line(106) = {106,101};
+Line(107) = {106,103};
+
+
+Line Loop(108) = {101, 102, -107, 106};
+Plane Surface(109) = {108};
+Line Loop(110) = {103, 104, 105, 107};
+Plane Surface(111) = {110};
+Physical Surface(112) = {109, 111};
+
+Transfinite Surface "*";
+Recombine Surface "*";
+Physical Surface(113) = {111, 109};
+
+Physical Line("XBlocked") = {103, 102};
+Physical Line("ImposedVelocity") = {105, 106};
+Physical Line("YBlocked") = {104, 101};
diff --git a/examples/python/custom-material-dynamics/material.dat b/examples/python/custom-material-dynamics/material.dat
new file mode 100644
index 000000000..cbb2e7b68
--- /dev/null
+++ b/examples/python/custom-material-dynamics/material.dat
@@ -0,0 +1,4 @@
+material local_elastic [
+	 name = steel
+	 rho = 1      # density
+]
diff --git a/examples/python/custom-material-dynamics/newmark.py b/examples/python/custom-material-dynamics/newmark.py
new file mode 100644
index 000000000..a7efcf261
--- /dev/null
+++ b/examples/python/custom-material-dynamics/newmark.py
@@ -0,0 +1,173 @@
+#!/usr/bin/python
+################################################################
+import akantu
+import numpy as np
+import os,subprocess
+################################################################
+class FixedValue:
+
+    def __init__(self,value,axis):
+        self.value = value
+        if axis == 'x': axis = 0
+        if axis == 'y': axis = 1
+        self.axis = axis
+        
+    def operator(self,node,flags,disp,coord):
+        # sets the displacement to the desired value in the desired axis
+        disp[self.axis] = self.value
+        # sets the blocked dofs vector to true in the desired axis
+        flags[self.axis] = True
+
+################################################################
+
+        
+class LocalElastic:
+
+    def __init__(self):
+
+        ## young modulus
+        self.E = 1
+        ## Poisson coefficient
+        self.nu = 0.3
+        ## density
+        self.rho = 1
+        ## First Lame coefficient
+        self._lambda = self.nu * self.E / ((1 + self.nu) * (1 - 2*self.nu))
+        ## Second Lame coefficient (shear modulus)
+        self.mu = self.E / (2 * (1 + self.nu));
+
+    ## declares all the internals
+    def registerInternals(self):
+        return []
+
+    ## declares all the parameters that could be parsed
+    def registerParam(self):
+        return []
+
+        
+    ## declares all the parameters that are needed
+    def getPushWaveSpeed(self):
+        return np.sqrt((self._lambda + 2*self.mu)/self.rho);
+
+    ## constitutive law for a given quadrature point
+    def computeStress(self,grad_u,sigma,internals):
+        lbda = 1.
+        mu = 1.
+        
+        trace = grad_u.trace();
+        sigma[:,:] = lbda*trace*np.eye(2) + mu * (grad_u + grad_u.T)
+
+
+################################################################    
+def main():
+
+    spatial_dimension = 2
+    Lbar = 10.
+    akantu.initialize('material.dat')
+    
+    mesh_file = 'bar.msh'
+    max_steps = 250
+    time_step = 1e-3
+
+    #if mesh was not created the calls gmsh to generate it
+    if not os.path.isfile(mesh_file):
+        import subprocess
+        ret = subprocess.call('gmsh -2 bar.geo bar.msh',shell=True)
+        if not ret == 0:
+            raise Exception('execution of GMSH failed: do you have it installed ?')
+
+    
+    ################################################################
+    ## Initialization
+    ################################################################
+    mesh = akantu.Mesh(spatial_dimension)
+    mesh.read(mesh_file)
+
+    mesh.createGroupsFromStringMeshData("physical_names")
+    model = akantu.SolidMechanicsModel(mesh)
+
+    model.initFull(akantu.SolidMechanicsModelOptions(akantu._explicit_lumped_mass,True))
+    mat = LocalElastic()
+    model.registerNewPythonMaterial(mat,"local_elastic")
+    model.initMaterials()
+
+    
+    model.setBaseName("waves")
+    model.addDumpFieldVector("displacement")
+    model.addDumpFieldVector("acceleration")
+    model.addDumpFieldVector("velocity")
+    model.addDumpField("blocked_dofs")
+
+    ################################################################
+    ## Boundary conditions
+    ################################################################
+    residual     = model.getResidual()
+    mass         = model.getMass()
+
+    displacement = model.getDisplacement()
+    acceleration = model.getAcceleration()
+    velocity     = model.getVelocity()
+
+    blocked_dofs = model.getBlockedDOFs()
+
+    ################################################################
+    ## boundary conditions
+    ################################################################
+    
+    model.applyDirichletBC(FixedValue(0,'x'), "XBlocked")
+    model.applyDirichletBC(FixedValue(0,'y'), "YBlocked")    
+
+    ################################################################
+    ## initial conditions
+    ################################################################
+
+    nb_nodes = mesh.getNbNodes()
+    position = mesh.getNodes()
+
+    pulse_width = 1
+    A = 0.01
+    for i in range(0,nb_nodes):
+        # Sinus * Gaussian
+        x = position[i, 0] - 5.
+        L = pulse_width
+        k = 0.1 * 2 * np.pi * 3 / L
+        displacement[i, 0] = A * np.sin(k * x) * np.exp(-(k * x) * (k * x) / (L * L))
+
+    ################################################################
+    ## timestep value computation
+    ################################################################
+    time_factor = 0.8
+    stable_time_step = model.getStableTimeStep() * time_factor
+
+    print "Stable Time Step = {0}".format(stable_time_step)
+    print "Required Time Step = {0}".format(time_step)
+
+    time_step = stable_time_step * time_factor
+
+    model.setTimeStep(time_step)
+
+    ################################################################
+    ## loop for evolution of motion dynamics
+    ################################################################
+    model.updateResidual()
+
+    epot = model.getEnergy('potential')
+    ekin = model.getEnergy('kinetic')
+
+    print "step,step * time_step,epot,ekin,epot + ekin"
+    for step in range(0,max_steps+1):
+        
+        model.dump()
+        ## output energy calculation to screen
+        print "{0},{1},{2},{3},{4}".format(step,step * time_step,
+                                           epot,ekin,
+                                           (epot + ekin))
+
+        model.solveStep()
+
+    akantu.finalize()
+    return 
+
+################################################################
+if __name__ == "__main__":
+    main()
diff --git a/examples/python/plate-hole/material.dat b/examples/python/plate-hole/material.dat
new file mode 100644
index 000000000..8f20aa1e1
--- /dev/null
+++ b/examples/python/plate-hole/material.dat
@@ -0,0 +1,6 @@
+material elastic [
+         name = steel
+	 rho = 7800     # density
+	 E   = 2.1e11   # young's modulus
+	 nu  = 0.3      # poisson's ratio
+]
diff --git a/examples/python/plate-hole/plate.geo b/examples/python/plate-hole/plate.geo
new file mode 100644
index 000000000..635cd9d48
--- /dev/null
+++ b/examples/python/plate-hole/plate.geo
@@ -0,0 +1,22 @@
+h1 = 0.1;
+h2 = h1;
+Point(1) = {0, 0, 0, h1};
+Point(2) = {4, 0, 0, h2};
+Point(3) = {4, 4, 0, h2};
+Point(4) = {0, 4, 0, h2};
+Point(5) = {1, 0, 0, h1};
+Point(6) = {0, 1, 0, h1};
+Circle(1) = {5, 1, 6};
+Line(2) = {6, 4};
+Line(3) = {4, 3};
+Line(4) = {3, 2};
+Line(5) = {2, 5};
+Line Loop(6) = {-1, -2, -3, -4, -5};
+Plane Surface(7) = {6};
+Physical Surface(8) = {7};
+Physical Line("XBlocked") = {2};
+Physical Line("YBlocked") = {5};
+Physical Line("Traction") = {3};
+
+Physical Point("XBlocked") = {5};
+Physical Point("YBlocked") = {5};
diff --git a/examples/python/plate-hole/plate.py b/examples/python/plate-hole/plate.py
new file mode 100644
index 000000000..d02f6e181
--- /dev/null
+++ b/examples/python/plate-hole/plate.py
@@ -0,0 +1,107 @@
+import akantu
+import numpy as np
+
+################################################################
+## Dirichlet Boudary condition functor: fix the displacement
+################################################################
+class FixedValue:
+
+    def __init__(self,value,axis):
+        self.value = value
+        if axis == 'x': axis = 0
+        if axis == 'y': axis = 1
+        self.axis = axis
+        
+    def operator(self,node,flags,disp,coord):
+        # sets the displacement to the desired value in the desired axis
+        disp[self.axis] = self.value
+        # sets the blocked dofs vector to true in the desired axis
+        flags[self.axis] = True
+
+################################################################
+## Neumann Boudary condition functor: impose a traction
+################################################################
+class FromTraction:
+
+    def __init__(self,traction):
+        self.traction = traction
+
+    def operator(self,quad_point,force,coord,normals):
+        # sets the force to the desired value in the desired axis
+        force[:] = self.traction
+
+################################################################
+
+def solve(material_file,mesh_file,traction):
+    akantu.initialize(material_file)
+    spatial_dimension = 2
+    
+    ################################################################
+    ##Initialization                                                          
+    ################################################################
+    mesh = akantu.Mesh(spatial_dimension)
+    mesh.read(mesh_file)
+    akantu.MeshUtils.purifyMesh(mesh)
+    mesh.createGroupsFromStringMeshData("physical_names")
+
+    model = akantu.SolidMechanicsModel(mesh)
+    model.initFull(akantu.SolidMechanicsModelOptions(akantu._static))
+    model.assembleStiffnessMatrix()
+    model.updateResidual()
+    
+    model.setBaseName("plate")
+    model.addDumpFieldVector("displacement")
+    model.addDumpFieldVector("force")
+    model.addDumpField("boundary")
+    model.addDumpField("strain")
+    model.addDumpField("stress")
+    model.addDumpField("blocked_dofs")
+    
+    ################################################################
+    ##Boundary conditions                                                      
+    ################################################################
+    displacement = model.getDisplacement()
+    blocked_dofs = model.getBlockedDOFs()
+    
+    model.applyDirichletBC(FixedValue(0.0, 'x'), "XBlocked")
+    model.applyDirichletBC(FixedValue(0.0, 'y'), "YBlocked")
+
+    trac = np.zeros(spatial_dimension)
+    trac[1] = traction
+
+    print "Solve for traction ",traction
+
+    model.getForce()[:] = 0
+    model.applyNeumannBC(FromTraction(trac), "Traction")
+
+    model.solveStaticDisplacement(1e-10,2);
+
+    model.dump()
+    akantu.finalize()
+
+################################################################
+# main 
+################################################################
+
+def main():
+
+    import os
+    mesh_file = 'plate.msh'
+    #if mesh was not created the calls gmsh to generate it
+    if not os.path.isfile(mesh_file):
+        import subprocess
+        ret = subprocess.call('gmsh -2 plate.geo {0}'.format(mesh_file),shell=True)
+        if not ret == 0:
+            raise Exception('execution of GMSH failed: do you have it installed ?')
+
+        
+    material_file = 'material.dat'
+    spatial_dimension = 2
+
+    traction = 1.
+    solve(material_file,mesh_file,traction)
+
+################################################################
+if __name__ == "__main__":
+    main() 
+