diff --git a/dynamical_system.py b/dynamical_system.py
index 3b60f8f..6291ac8 100644
--- a/dynamical_system.py
+++ b/dynamical_system.py
@@ -312,7 +312,7 @@ def propagate(self, nu1, nu2, x1):
                 x2_bar = None
                 if self.mu != 0 and (self.Li == 1 or self.Li == 2 or self.Li == 3):
                     if self.ecc == 0.:
-                        phi = phi_harmo(puls_oop_LP(self.x_L_normalized, self.mu) * (nu2 - nu1))
+                        phi = phi_harmo(nu2 - nu1, puls_oop_LP(self.x_L_normalized, self.mu))
                         x2_bar = phi . dot(x1_bar)
                     else:  # elliptical case
                         print('PROPAGATE: 3-body elliptical out-of-plane near L1, 2 and 3 not coded yet')
@@ -356,13 +356,14 @@ def compute_rhs(self, BC):
         if BC.half_dim == 1:
             if (self.mu != 0.) and (self.Li == 1 or self.Li == 2 or self.Li == 3):
                 if self.ecc == 0.:
-                    u += phi_harmo(-puls_oop_LP(self.x_L_normalized, self.mu) * BC.nuf) . dot(x2)
-                    u -= phi_harmo(-puls_oop_LP(self.x_L_normalized, self.mu) * BC.nu0) . dot(x1)
+                    pulsation = puls_oop_LP(self.x_L_normalized, self.mu)
+                    u += phi_harmo(-BC.nuf, pulsation). dot(x2)
+                    u -= phi_harmo(-BC.nu0, pulsation). dot(x1)
                 else:
                     print('compute_rhs: 3-body elliptical out-of-plane near L1, 2 and 3 not coded yet')
             else:  # out-of-plane elliptical 2-body problem or 3-body near L4 and 5
-                u += phi_harmo(-BC.nuf).dot(x2)
-                u -= phi_harmo(-BC.nu0).dot(x1)
+                u += phi_harmo(-BC.nuf, 1.0).dot(x2)
+                u -= phi_harmo(-BC.nu0, 1.0).dot(x1)
             u[0] *= multiplier
             u[1] *= multiplier
         elif BC.half_dim == 2:
diff --git a/indirect_solver.py b/indirect_solver.py
index c60e2ac..9680c7f 100644
--- a/indirect_solver.py
+++ b/indirect_solver.py
@@ -116,10 +116,8 @@ def _circular_oop_L123(self, BC):
         puls = orbital_mechanics.puls_oop_LP(self.dyn.x_L_normalized, self.dyn.mu)
 
         # scaling inputs
-        x0_rescaled = numpy.array((BC.x0[0], BC.x0[1] / puls))
-        xf_rescaled = numpy.array((BC.xf[0], BC.xf[1] / puls))
-        BC_rescaled = utils.BoundaryConditions(BC.nu0 * puls, BC.nuf * puls, x0_rescaled, xf_rescaled)
-        dyn_rescaled = dynamical_system.DynamicalSystem(0., 0., self.dyn.period * puls, self.dyn.Li)
+        BC_rescaled = utils.BoundaryConditions(BC.nu0 * puls, BC.nuf * puls, BC.x0, BC.xf)
+        dyn_rescaled = dynamical_system.DynamicalSystem(0., 0., self.dyn.period / puls, self.dyn.sma)
         z_rescaled = dyn_rescaled.compute_rhs(BC_rescaled)
 
         # 'classic' call to numerical solver
@@ -128,6 +126,8 @@ def _circular_oop_L123(self, BC):
         # un-scaling
         for k in range(0, len(nus)):
             nus[k] /= puls
-            DVs[k] *= puls
+        factor = linalg.norm(self.dyn.compute_rhs(BC)) / linalg.norm(dyn_rescaled.compute_rhs(BC_rescaled))
+        for k in range(0, len(lamb)):
+            lamb[k] /= factor
 
         return utils.ControlLaw(BC.half_dim, nus, DVs, lamb)
diff --git a/moments.py b/moments.py
index 32fe4ed..431467b 100644
--- a/moments.py
+++ b/moments.py
@@ -25,7 +25,7 @@ def Y_oop(e, nu):
     """
 
     rho = rho_func(e, nu)
-    phi = phi_harmo(-nu)
+    phi = phi_harmo(-nu, 1.0)
     Y = numpy.zeros((2, 1))
     Y[0, 0] = phi[0, 1] / rho
     Y[1, 0] = phi[1, 1] / rho
@@ -47,7 +47,7 @@ def Y_oop_LP123(nu, x, mu):
 
     """
 
-    phi = phi_harmo(-puls_oop_LP(x, mu) * nu)
+    phi = phi_harmo(-nu, puls_oop_LP(x, mu))
     Y = numpy.zeros((2, 1))
     Y[0, 0] = phi[0, 1]
     Y[1, 0] = phi[1, 1]
diff --git a/orbital_mechanics.py b/orbital_mechanics.py
index 32c311c..7ebf938 100644
--- a/orbital_mechanics.py
+++ b/orbital_mechanics.py
@@ -321,11 +321,12 @@ def rho_func(e, nu):
     return 1.0 + e * math.cos(nu)
 
 
-def phi_harmo(nu):
+def phi_harmo(nu, pulsation):
     """Function returning the transition matrix of the harmonic oscillator.
 
             Args:
                 nu (float): true anomaly.
+                pulsation (float): pulsation.
 
             Returns:
                 phi (numpy.array): transition matrix of harmonic oscillator.
@@ -333,10 +334,12 @@ def phi_harmo(nu):
     """
 
     phi = numpy.zeros((2, 2))
-    phi[0, 0] = math.cos(nu)
-    phi[0, 1] = math.sin(nu)
-    phi[1, 0] = -phi[0, 1]
-    phi[1, 1] = phi[0, 0]
+    c = math.cos(pulsation * nu)
+    s = math.sin(pulsation * nu)
+    phi[0, 0] = c
+    phi[0, 1] = s / pulsation
+    phi[1, 0] = -s * pulsation
+    phi[1, 1] = c
 
     return phi
 
@@ -358,7 +361,7 @@ def transition_oop(x1_bar, nu1, nu2):
     if len(x1_bar) != 2:
         print('TRANSITION_OOP: out-of-plane initial conditions need to be two-dimensional')
 
-    return phi_harmo(nu2 - nu1) . dot(x1_bar)
+    return phi_harmo(nu2 - nu1, 1.0). dot(x1_bar)
 
 
 def exp_HCW(nu):