--- a
+++ b/inst/example/Elasticity/HyperElasticity.m
@@ -0,0 +1,57 @@
+pkg load fem-fenics msh
+problem = 'HyperElasticity';
+import_ufl_Problem (problem);
+
+Rotation = @(x,y,z) ...
+   [0; ...
+   0.5*(0.5 + (y - 0.5)*cos(pi/3) - (z-0.5)*sin(pi/3) - y);...
+   0.5*(0.5 + (y - 0.5)*sin(pi/3) + (z-0.5)*cos(pi/3) - z)];
+
+#Create mesh and define function space
+x = y = z = linspace (0, 1, 17);
+mshd = Mesh (msh3m_structured_mesh (x, y, z, 1, 1:6));
+V  = FunctionSpace (problem, mshd);
+
+# Create Dirichlet boundary conditions
+bcl = DirichletBC (V, @(x,y,z) [0; 0; 0], 1);
+bcr = DirichletBC (V, Rotation, 2);
+bcs = {bcl, bcr};
+ 
+# Define source and boundary traction functions
+B = Constant ('B', [0.0; -0.5; 0.0]);
+T = Constant ('T', [0.1; 0.0; 0.0]);
+
+# Set material parameters
+E = 10.0;
+nu = 0.3;
+mu = Constant ('mu', E./(2*(1+nu)));
+lmbda = Constant ('lmbda', E*nu./((1+nu)*(1-2*nu)));
+u = Expression ('u', @(x,y,z) [0; 0; 0]);
+
+# Create (linear) form defining (nonlinear) variational problem
+L = ResidualForm (problem, V, mu, lmbda, B, T, u);
+
+# Solve nonlinear variational problem F(u; v) = 0
+
+# Create function for the resolution of the NL problem
+  function [y, jac] = f (problem, xx, V, bc1, bc2, B, T, mu, lmbda)
+    u = Function ('u', V, xx);
+    a = JacobianForm (problem, V, mu, lmbda, u);
+    L = ResidualForm (problem, V, mu, lmbda, B, T, u);
+    if (nargout == 1)
+      [y, xx] = assemble (L, xx, bc1, bc2);
+    elseif (nargout == 2)
+      [jac, y, xx] = assemble_system (a, L, xx, bc1, bc2);
+    endif
+  endfunction
+
+u0 = assemble (L, bcs{:});
+fs = @(xx) f (problem, xx, V, bcl, bcr, B, T, mu, lmbda);
+[x, fval, info] = fsolve (fs, u0, optimset ("jacobian", "on"));
+func = Function ('u', V, x);
+ 
+# Save solution in VTK format
+save (func, 'displacement');
+
+# Plot solution
+plot (func);