diff --git a/bin/args.py b/bin/args.py
index ae606d86540a344b1e1b7d7bf037a1d2dd07f294..013b0db39febdea04f7972627489d51d1711c27a 100644
--- a/bin/args.py
+++ b/bin/args.py
@@ -10,7 +10,7 @@ See the LICENSE.txt and README.md for more license and copyright information.
 import argparse
 
 def parse():
-    """Parses the python script arguments
+    """Parses python script arguments
 
     Returns an argparse
     """
@@ -69,8 +69,7 @@ desc = ("a non-linear eigenvalue solver " +
 
 
 class MyFormatter(argparse.HelpFormatter):
-    """ Custom Formatter
-    """
+    """ Custom Formatter"""
 
     def __init__(self, prog):
         super(MyFormatter, self).__init__(prog,
diff --git a/bin/beyn.py b/bin/beyn.py
index 177d33c2b235a87856bd4edd0c2af4095904c6ca..3129d80548ab92f8f308fa03cee51d6ee1e4ba3f 100644
--- a/bin/beyn.py
+++ b/bin/beyn.py
@@ -44,9 +44,6 @@ class Beyn:
         self.__rankTol  = rankTol
         self.__verbose  = verbose
 
-        # Display the parameters
-        if(verbose): self.__display()
-
     def simple(self):
         """Solves the eigenvalue problem using Beyn's algorithm (simple version)
 
@@ -131,6 +128,32 @@ class Beyn:
         if(verbose): print "Done!"
         return myLambda, Q, norm
 
+    def quite(self, beQuite):
+        """Sets this Beyn solver verbosity
+
+        If beQuite is True, this solver becomes quite; it is verbose otherwise
+        """
+        self.__verbose = not beQuite
+
+    def display(self):
+        """Prints this Beyn solver parameters"""
+        print "Beyn's contour integral method (simple)"
+        print "---------------------------------------"
+        print " # Nodes used for the trapezoidal rule:"+ " " +str(self.__nNode)
+        print " # Maximum number of iterations:       "+ " " +str(self.__maxIt)
+        print " # Initial size of col(A0):            "+ " " +str(self.__lStart)
+        print " # Step size for increasing col(A0):   "+ " " +str(self.__lStep)
+        print " # Rank test relative tolerance:       ",
+        print format(self.__rankTol, ".2e")
+        print "---------------------------------------"
+        print " # Cirular path origin:                ",
+        print "(" + format(np.real(self.__origin), "+.2e") + ")",
+        print "+",
+        print "(" + format(np.imag(self.__origin), "+.2e") + ")j"
+        print " # Cirular path radius:                ",
+        print format(self.__radius, "+.2e")
+        print "---------------------------------------"
+
     def __rank(self, A):
         """Returns the rank of a given matrix and its singular values"""
         S   = np.linalg.svd(A, full_matrices=False, compute_uv=0)
@@ -195,22 +218,3 @@ class Beyn:
     def __randomMatrix(n, m):
         """Returns a random complex matrix of the given size (n x m)"""
         return np.matlib.rand(n, m) + np.matlib.rand(n, m) * 1j
-
-    def __display(self):
-        """Prints this Beyn parameters"""
-        print "Beyn's contour integral method (simple)"
-        print "---------------------------------------"
-        print " # Nodes used for the trapezoidal rule:"+ " " +str(self.__nNode)
-        print " # Maximum number of iterations:       "+ " " +str(self.__maxIt)
-        print " # Initial size of col(A0):            "+ " " +str(self.__lStart)
-        print " # Step size for increasing col(A0):   "+ " " +str(self.__lStep)
-        print " # Rank test relative tolerance:       ",
-        print format(self.__rankTol, ".2e")
-        print "---------------------------------------"
-        print " # Cirular path origin:                ",
-        print "(" + format(np.real(self.__origin), "+.2e") + ")",
-        print "+",
-        print "(" + format(np.imag(self.__origin), "+.2e") + ")j"
-        print " # Cirular path radius:                ",
-        print format(self.__radius, "+.2e")
-        print "---------------------------------------"
diff --git a/bin/cim.py b/bin/cim.py
index 4973fc073ef76c29cf37f66acca47aa5e706b755..bcc4acc51b7f9d3cf77b89ef08d9cfdda3e268da 100755
--- a/bin/cim.py
+++ b/bin/cim.py
@@ -8,11 +8,12 @@ Technische Universitaet Darmstadt.
 See the LICENSE.txt and README.md for more license and copyright information.
 """
 
-import solver as Solver
-import beyn   as Beyn
-import numpy  as np
-import args   as args
-import time   as timer
+from solver       import GetDPWave
+from beyn         import Beyn
+from fixed_point  import FixedPoint
+import args       as args
+import numpy      as np
+import time       as timer
 
 # Start timer
 tStart = timer.time()
@@ -28,18 +29,25 @@ for i in arg.setnumber:
         setnumber.append(j)
 
 # Initialise GetDP
-operator = Solver.GetDPWave(arg.pro, arg.mesh, arg.resolution, setnumber)
+operator = GetDPWave(arg.pro, arg.mesh, arg.resolution, setnumber)
 
 # Initialise Beyn's solver
-beyn = Beyn.Beyn(operator, arg.origin, arg.radius,
-                 arg.nodes, arg.maxIt, arg.lStart, arg.lStep, arg.rankTol,
-                 not arg.quiet)
+beyn = Beyn(operator, arg.origin, arg.radius,
+            arg.nodes, arg.maxIt, arg.lStart, arg.lStep, arg.rankTol,
+            not arg.quiet)
+
+# Display parameters
+if(not arg.quiet):
+    beyn.display()
 
 # Compute
 if(not arg.non_holomorphic):
     l, V, r = beyn.simple()
 else:
-    raise ValueError('Non-holomorphic case not implemented... yet')
+    fixed = FixedPoint(beyn) # This variable ('fixed') is meant to be called
+    fixed.solve()            # elsewhere (i.e. Python[] call in GetDP) thanks to
+                             # the python interpreter.
+                             # Please don't change its name!
 
 # Display the computed eigenvalues
 if(not arg.quiet):
diff --git a/bin/fixed_point.py b/bin/fixed_point.py
new file mode 100644
index 0000000000000000000000000000000000000000..3df7718053c903808be6859cfcca5d73f4f494be
--- /dev/null
+++ b/bin/fixed_point.py
@@ -0,0 +1,64 @@
+"""
+cim.py, a non-linear eigenvalue solver.
+Copyright (C) 2017 N. Marsic, F. Wolf, S. Schoeps and H. De Gersem,
+Institut fuer Theorie Elektromagnetischer Felder (TEMF),
+Technische Universitaet Darmstadt.
+
+See the LICENSE.txt and README.md for more license and copyright information.
+"""
+
+import beyn as Beyn
+
+"""
+At some point I will need:
+*/ a cluster identification
+*/ an error with respect to the *NON*-holomorphic problem
+*/ a way to switch between holomorphic and non-holomorphic runs
+*/ a corrector?
+"""
+
+class FixedPoint:
+    """A fixed-point solver for handling non-holomorphic operator
+
+    The fixed-point approach follows the strategy bellow:
+    TODO
+
+    if cluster = [], use corrector initial guess
+    if needTrueResidual = True, corrector becomes non-holomorphic
+    """
+
+    def __init__(self, solver):
+        """Instantiates a new FixedPoint non-holomorphic solver
+
+        Keyword arguments:
+        solver -- the solver used to compute holomorphic solutions
+        """
+        # Save
+        self.__solver = solver
+
+        # State
+        self.__cluster     = []
+        self.__holomorphic = True
+
+    def solve(self):
+        """Solves the eigenvalue problem by a fixed-point iteration
+
+        It returns the computed eigenvalues, eigenvectors, TODO
+        """
+        # (Re)initialize
+        self.__cluster     = []
+        self.__holomorphic = True
+
+        # First run with corrector initial guess
+        l, V, r = self.__solver.simple()
+
+        print l
+
+    def cluster(self):
+        """Returns the eigenvalues (by cluster) found by the last iteration"""
+        return self.__cluster
+
+    def needNonHolomorphic(self):
+        """This method returns True if this FixedPoint needs the true,
+        non-holomorphic, operator and not its holomorphic approximation"""
+        return self.__trueRes
diff --git a/bin/interpolator.py b/bin/interpolator.py
new file mode 100644
index 0000000000000000000000000000000000000000..61d63aa97e11e668f3acbd1979babefe8345d10a
--- /dev/null
+++ b/bin/interpolator.py
@@ -0,0 +1,33 @@
+"""
+cim.py, a non-linear eigenvalue solver.
+Copyright (C) 2017 N. Marsic, F. Wolf, S. Schoeps and H. De Gersem,
+Institut fuer Theorie Elektromagnetischer Felder (TEMF),
+Technische Universitaet Darmstadt.
+
+See the LICENSE.txt and README.md for more license and copyright information.
+"""
+
+from numpy import poly1d
+from scipy import interpolate
+from abc   import ABCMeta, abstractmethod
+
+class Interpolator:
+    """A generic class for interpolation"""
+    @abstractmethod
+    def __init__(self):
+        pass
+
+    def at(self, x):
+        """Returns the value of this Interpolation at a given point"""
+        return self._p(x) # We use poly1d from scipy
+
+    __metaclass__ = ABCMeta
+
+
+class Lagrange(Interpolator):
+    """A class for Lagrange interpolation"""
+    def __init__(self, xdata, ydata):
+        """Instantiates a new Interpolator, such that:
+        p(xdata) = ydata for all xdata, where p is build with a Lagrange basis
+        """
+        self._p = interpolate.lagrange(xdata, ydata)
diff --git a/bin/solver.py b/bin/solver.py
index f52eec5a439da3c88de0c214b59d5bc0421165ed..c6e7254f62d5d674cae4828410eb4f60cc4bde6e 100644
--- a/bin/solver.py
+++ b/bin/solver.py
@@ -62,7 +62,7 @@ class GetDPWave:
         return self.solution(0).shape[0]
 
     def apply(self, x, w):
-        """Applies x to the operator with a pulsation of w
+        """Applies x to the operator with an angular frequency of w
 
         This method updates self.solution()
         """
diff --git a/example/maxwell_sibc_lagrange/Makefile b/example/maxwell_sibc_lagrange/Makefile
index b7a7c71157fba5c6f0b2ad2f0e959e49ae93dbfb..76d2e4439a5f22f43e3f9f069a761a73406b37b3 100644
--- a/example/maxwell_sibc_lagrange/Makefile
+++ b/example/maxwell_sibc_lagrange/Makefile
@@ -6,12 +6,12 @@ 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 1e10+1e9j 4e9 -nodes 30 -lStart 20 -lStep 30 -setnumber isSC 0 -non_holomorphic
 
-#	python ${BIN}/cim.py sphere.pro sphere.msh Solve 7.4e9+7.3e8j 1e9 -nodes 10 -lStart 4 -setnumber isSC 0
+	python ${BIN}/cim.py sphere.pro sphere.msh Solve 7.4e9+7.3e8j 1e9 -nodes 10 -lStart 4 -setnumber isSC 0 -non_holomorphic
 #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
+#	python ${BIN}/cim.py sphere.pro sphere.msh Solve 1.0e10+1.1e9j 1e9 -nodes 10 -lStart 6 -setnumber isSC 0 -non_holomorphic
 #1.050336615299e+10 + 1.186115651023e+09j
 
 clean:
diff --git a/example/maxwell_sibc_lagrange/cimResolution.pro b/example/maxwell_sibc_lagrange/cimResolution.pro
index 99551870563d2775702a773028f245864013abaa..033844d44d1dd5796d1d810e84c61a4b0c116706 100644
--- a/example/maxwell_sibc_lagrange/cimResolution.pro
+++ b/example/maxwell_sibc_lagrange/cimResolution.pro
@@ -7,6 +7,10 @@ Resolution{
   { Name Solve;
     System{ { Name A; NameOfFormulation FormulationCIM; Type ComplexValue; } }
     Operation{
+      If(!doInit)
+        Evaluate[Python[]{"real_corrector_init.py"}];
+      EndIf
+
       Generate[A];
 
       If(doInit)
diff --git a/example/maxwell_sibc_lagrange/real_corrector_handl.py b/example/maxwell_sibc_lagrange/real_corrector_handl.py
new file mode 100644
index 0000000000000000000000000000000000000000..359e7b42aa3474eb0d7bf74b304ab43fd2d16357
--- /dev/null
+++ b/example/maxwell_sibc_lagrange/real_corrector_handl.py
@@ -0,0 +1,20 @@
+"""Real-part operator corrector
+
+Approximates the real-part operator by an holomorphic function, constructed with
+a Lagrange polynomial defined at some correction points. Please make sure that
+GetDP calls real_corrector_init.py before calling this file.
+
+The naming convention below must be followed:
+-- a FixedPoint object named 'fixed' is available and prepared by cim.py
+-- a Lagrange Interpolator object named 'lagrange' coming from initialization
+-- the variable 'output' is read back by GetDP as the Pyhton[] call output
+"""
+
+# Input angular frequency
+omega = input[0]
+
+# Compute holomorphic approximation of the real part operator
+if(not fixed.cluster()):
+    output = omega                            # If no cluster: Re(omega) ~ omega
+else:
+    output = (omega + lagrange.at(omega)) / 2 # Else: correct
diff --git a/example/maxwell_sibc_lagrange/real_corrector_init.py b/example/maxwell_sibc_lagrange/real_corrector_init.py
new file mode 100644
index 0000000000000000000000000000000000000000..942aa3b9bd16b11c4fc6a22250d18ece71eb4cb1
--- /dev/null
+++ b/example/maxwell_sibc_lagrange/real_corrector_init.py
@@ -0,0 +1,13 @@
+"""Real-part operator corrector: initialization
+
+See: real_corrector_handl.py for more information
+"""
+
+from interpolator import Lagrange
+
+# Previous eigenvalue clusters
+cluster = fixed.cluster()
+
+# Initialize interpolator
+if(cluster):
+    lagrange = Lagrange(cluster, np.conj(cluster).tolist())
diff --git a/example/maxwell_sibc_lagrange/sphere.pro b/example/maxwell_sibc_lagrange/sphere.pro
index a044bb8afe20bd7c672ef46ed3c4667e8ec2af05..a7c28b0d7c3e16a58be5d9aff7929aec2e1c29ed 100644
--- a/example/maxwell_sibc_lagrange/sphere.pro
+++ b/example/maxwell_sibc_lagrange/sphere.pro
@@ -16,27 +16,24 @@ Function{
   Mu0  = 4*Pi*1e-7;      // [H/m]
   Eps0 = 1/(Mu0 * C0^2); // [F/m]
 
+  // Correction for real-part operator //
+  If(doInit)
+    PyPulsRe[] = 1;
+  EndIf
+  If(!doInit)
+    PyPulsRe[] = Python[Puls[]]{"real_corrector_handl.py"};
+  EndIf
+
   // Conductivity (coming from sphere.dat) //
-  If(isSC == 1)                                            // Superconductor
-    Sigma[] = Sigma - J[] / (Mu0 * Lambda^2 * Re[Puls[]]); // [S/m]
+  If(isSC == 1)                                             // Superconductor
+    Sigma[] = Sigma - J[] / (Mu0 * Lambda^2 * PyPulsRe[]]); // [S/m]
   EndIf
-  If(isSC == 0)                                            // Normal conductor
-    Sigma[] = Sigma;                                       // [S/m]
+  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
+  SL[] = Sqrt[(J[] * PyPulsRe[] * Mu0) / Sigma[]]; // H-Formulation
   SIbc[] = -1.0 / SL[];                            // E-Formulation
 }