diff --git a/example/maxwell_sibc_lagrange/Makefile b/example/maxwell_sibc_lagrange/Makefile
new file mode 100644
index 0000000000000000000000000000000000000000..b7a7c71157fba5c6f0b2ad2f0e959e49ae93dbfb
--- /dev/null
+++ b/example/maxwell_sibc_lagrange/Makefile
@@ -0,0 +1,21 @@
+BIN=../../bin
+
+sphere: init cim
+
+init:
+ gmsh sphere.geo -3
+
+cim:
+ python ${BIN}/cim.py sphere.pro sphere.msh Solve 1e10+1e9j 4e9 -nodes 20 -lStart 20 -setnumber isSC 0
+
+# python ${BIN}/cim.py sphere.pro sphere.msh Solve 7.4e9+7.3e8j 1e9 -nodes 10 -lStart 4 -setnumber isSC 0
+#7.474355560164e+09 + 7.742536715699e+08j
+
+# python ${BIN}/cim.py sphere.pro sphere.msh Solve 1.0e10+1.1e9j 1e9 -nodes 10 -lStart 6 -setnumber isSC 0
+#1.050336615299e+10 + 1.186115651023e+09j
+
+clean:
+ rm -f *.pos
+ rm -f *.msh
+ rm -f *.pre
+ rm -f *.res
diff --git a/example/maxwell_sibc_lagrange/cimParameters.pro b/example/maxwell_sibc_lagrange/cimParameters.pro
new file mode 100644
index 0000000000000000000000000000000000000000..80b91c50abbb3e88012b971fd11ed04ad33c876a
--- /dev/null
+++ b/example/maxwell_sibc_lagrange/cimParameters.pro
@@ -0,0 +1,26 @@
+/////////////////////////////////////////////////////////////////
+// Parameters for contour integral method. //
+// The master file should inclue "cimResolution.pro" later on. //
+/////////////////////////////////////////////////////////////////
+Function{
+ // Physical data //
+ DefineConstant[angularFreqRe = 1, // [rad/s]
+ angularFreqIm = 0]; // [rad/s]
+
+ Puls[] = Complex[angularFreqRe,
+ angularFreqIm]; // [rad/m]
+
+ // Algebraic data //
+ DefineConstant[nRHS = 1]; // Number of RHS for this run
+ For i In {0:nRHS-1}
+ DefineConstant[x~{i}() = {}, // Solution
+ b~{i}() = {}]; // Right hand side
+ EndFor
+
+ // Control data //
+ DefineConstant[doInit = 0, // Should I initialize some stuff?
+ doSolve = 0, // Should I solve Ax = b?
+ doPostpro = 0, // Should I only create a view for x()?
+ doApply = 0, // Should I only apply x(): x <- Ax?
+ fileName = "eig.pos"]; // Postpro file name
+}
diff --git a/example/maxwell_sibc_lagrange/cimResolution.pro b/example/maxwell_sibc_lagrange/cimResolution.pro
new file mode 100644
index 0000000000000000000000000000000000000000..99551870563d2775702a773028f245864013abaa
--- /dev/null
+++ b/example/maxwell_sibc_lagrange/cimResolution.pro
@@ -0,0 +1,42 @@
+/////////////////////////////////////////////////////////////////////
+// Resolution for contour integral method. //
+// Warning, the master file should include "cimParameters.pro", //
+// and should also define "FormulationCIM" the formulation to use. //
+/////////////////////////////////////////////////////////////////////
+Resolution{
+ { Name Solve;
+ System{ { Name A; NameOfFormulation FormulationCIM; Type ComplexValue; } }
+ Operation{
+ Generate[A];
+
+ If(doInit)
+ CopySolution[A, x~{0}()];
+ EndIf
+
+ If(doSolve)
+ // Full solve for first RHS
+ CopyRightHandSide[b~{0}(), A];
+ Solve[A];
+ CopySolution[A, x~{0}()];
+
+ // SolveAgain for other RHSs
+ For i In {1:nRHS-1}
+ CopyRightHandSide[b~{i}(), A];
+ SolveAgain[A];
+ CopySolution[A, x~{i}()];
+ EndFor
+ EndIf
+
+ If(doApply)
+ CopySolution[x~{0}(), A];
+ Apply[A];
+ CopySolution[A, x~{0}()];
+ EndIf
+
+ If(doPostpro)
+ CopySolution[x~{0}(), A];
+ PostOperation[Post];
+ EndIf
+ }
+ }
+}
diff --git a/example/maxwell_sibc_lagrange/sphere.dat b/example/maxwell_sibc_lagrange/sphere.dat
new file mode 100644
index 0000000000000000000000000000000000000000..db37ad8585b2209ae7880972c6fc1857206727ea
--- /dev/null
+++ b/example/maxwell_sibc_lagrange/sphere.dat
@@ -0,0 +1,22 @@
+// Geometry //
+DefineConstant[R = { 100e-3, Name "Input/00Geometry/00Radius [m]" }];
+
+// Mesh //
+DefineConstant[md = { 10, Name "Input/01Mesh/00Density [Radius^-1]" }];
+cl = R/md;
+
+// Material //
+DefineConstant[isSC = { 0,
+ Name "Input/02Material/00Is superconductor?",
+ GmshOption "Reset",
+ Choices {0 = "No", 1 = "Yes"} }];
+If(isSC == 0)
+DefineConstant[Sigma = { 1e00,
+ Name "Input/02Material/01Conductivity [S m^-1]" }];
+EndIf
+If(isSC == 1)
+DefineConstant[Sigma = { 1e9,
+ Name "Input/02Material/01Conductivity [S m^-1]" }];
+DefineConstant[Lambda = { 1e-8,
+ Name "Input/02Material/02London depth [m]" }];
+EndIf
diff --git a/example/maxwell_sibc_lagrange/sphere.geo b/example/maxwell_sibc_lagrange/sphere.geo
new file mode 100644
index 0000000000000000000000000000000000000000..5ad10a3b551aa48ff4a5da29ae48d2c965ec4cbc
--- /dev/null
+++ b/example/maxwell_sibc_lagrange/sphere.geo
@@ -0,0 +1,50 @@
+Include "sphere.dat";
+
+Point(0) = { 0, 0, 0, cl};
+
+Point(1) = {+R, 0, 0, cl};
+Point(2) = { 0, +R, 0, cl};
+Point(3) = {-R, 0, 0, cl};
+Point(4) = { 0, -R, 0, cl};
+Point(5) = { 0, 0, +R, cl};
+Point(6) = { 0, 0, -R, cl};
+
+Circle(01) = {1, 0, 2};
+Circle(02) = {2, 0, 3};
+Circle(03) = {3, 0, 4};
+Circle(04) = {4, 0, 1};
+
+Circle(05) = {5, 0, 1};
+Circle(06) = {1, 0, 6};
+Circle(07) = {6, 0, 3};
+Circle(08) = {3, 0, 5};
+
+Circle(09) = {2, 0, 6};
+Circle(10) = {6, 0, 4};
+Circle(11) = {4, 0, 5};
+Circle(12) = {5, 0, 2};
+
+Line Loop(1) = {1, -12, 5};
+Line Loop(2) = {1, 9, -6};
+Line Loop(3) = {9, 7, -2};
+Line Loop(4) = {2, 8, 12};
+Line Loop(5) = {11, 5, -4};
+Line Loop(6) = {4, 6, 10};
+Line Loop(7) = {10, -3, -7};
+Line Loop(8) = {3, 11, -8};
+
+Ruled Surface(1) = {1};
+Ruled Surface(2) = {2};
+Ruled Surface(3) = {3};
+Ruled Surface(4) = {4};
+Ruled Surface(5) = {5};
+Ruled Surface(6) = {6};
+Ruled Surface(7) = {7};
+Ruled Surface(8) = {8};
+
+Surface Loop(1) = {3, 2, 1, 4, 8, 7, 6, 5};
+Volume(1) = {1};
+
+Physical Volume(1) = {1};
+Physical Surface(2) = {1, 2, 3, 4, 5, 6, 7, 8};
+Physical Point(3) = {0};
diff --git a/example/maxwell_sibc_lagrange/sphere.pro b/example/maxwell_sibc_lagrange/sphere.pro
new file mode 100644
index 0000000000000000000000000000000000000000..a044bb8afe20bd7c672ef46ed3c4667e8ec2af05
--- /dev/null
+++ b/example/maxwell_sibc_lagrange/sphere.pro
@@ -0,0 +1,110 @@
+Include "sphere.dat";
+
+Group{
+ Omega = Region[1];
+ Bndry = Region[2];
+}
+
+Include "./cimParameters.pro"; // Functions for cim.py
+
+Function{
+ // Imaginary Unit //
+ J[] = Complex[0, 1];
+
+ // Physical Constants //
+ C0 = 299792458; // [m/s]
+ Mu0 = 4*Pi*1e-7; // [H/m]
+ Eps0 = 1/(Mu0 * C0^2); // [F/m]
+
+ // Conductivity (coming from sphere.dat) //
+ If(isSC == 1) // Superconductor
+ Sigma[] = Sigma - J[] / (Mu0 * Lambda^2 * Re[Puls[]]); // [S/m]
+ EndIf
+ If(isSC == 0) // Normal conductor
+ Sigma[] = Sigma; // [S/m]
+ EndIf
+
+ // Leontovich SIBC (H-Formulation) //
+ PulsG1[] = 7.474355560164e+09 + 7.742536715699e+08*J[];
+ PulsG2[] = 1.050336615299e+10 + 1.186115651023e+09*J[];
+
+ PulsBar1[] = Re[PulsG1[]] - J[]*Im[PulsG1[]];
+ PulsBar2[] = Re[PulsG2[]] - J[]*Im[PulsG2[]];
+
+ L1[] = (Puls[] - PulsG2[]) / (PulsG1[] - PulsG2[]);
+ L2[] = (Puls[] - PulsG1[]) / (PulsG2[] - PulsG1[]);
+
+ //PulsCor[] = PulsBar1[]*Puls[]/PulsG1[]; //Re[Puls[]] + J[]*Im[Puls[]];
+ PulsCor[] = PulsBar1[]*L1[] + PulsBar2[]*L2[];
+ SL[] = Sqrt[(J[] * (Puls[]+PulsCor[])/2 * Mu0) / Sigma[]]; // H-Formulation
+ SIbc[] = -1.0 / SL[]; // E-Formulation
+}
+
+Integration{
+ { Name I1;
+ Case{
+ { Type Gauss;
+ Case{
+ { GeoElement Tetrahedron; NumberOfPoints 4; } // O1 * O1
+ { GeoElement Triangle; NumberOfPoints 3; } // O1 * O1
+ { GeoElement Line; NumberOfPoints 2; } // O1 * O1
+ }
+ }
+ }
+ }
+}
+
+Jacobian{
+ { Name Jac;
+ Case{
+ { Region Omega; Jacobian Vol; }
+ { Region Bndry; Jacobian Sur; }
+ }
+ }
+}
+
+FunctionSpace{
+ { Name HCurl; Type Form1;
+ BasisFunction{
+ { Name sx; NameOfCoef ux; Function BF_Edge;
+ Support Region[{Omega, Bndry}]; Entity EdgesOf[All]; }
+ }
+ }
+}
+
+Formulation{
+ { Name FormulationCIM; Type FemEquation;
+ Quantity{
+ { Name e; Type Local; NameOfSpace HCurl; }
+ }
+ Equation{
+ // Maxwell
+ Galerkin{ [ 1/Mu0 * Dof{d e}, {d e}];
+ In Omega; Jacobian Jac; Integration I1; }
+ Galerkin{ [ -Puls[]^2 * Eps0 * Dof{ e}, { e}];
+ In Omega; Jacobian Jac; Integration I1; }
+
+ // IBC
+ Galerkin{ [-J[] * Puls[] * SIbc[] * Dof{ e}, { e}];
+ In Bndry; Jacobian Jac; Integration I1; }
+ }
+ }
+}
+
+Include "./cimResolution.pro"; // Resolution for cim.py
+
+PostProcessing{
+ { Name Post; NameOfFormulation FormulationCIM;
+ Quantity{
+ { Name E; Value{ Term{ [{e}]; In Omega; Jacobian Jac; } } }
+ }
+ }
+}
+
+PostOperation{
+ { Name Post; NameOfPostProcessing Post;
+ Operation{
+ Print[E, OnElementsOf Omega, File fileName];
+ }
+ }
+}