diff --git a/bioptim/optimization/solution.py b/bioptim/optimization/solution.py
index dcecebbd8..54d38c4ad 100644
--- a/bioptim/optimization/solution.py
+++ b/bioptim/optimization/solution.py
@@ -8,7 +8,7 @@
 from matplotlib import pyplot as plt
 
 from ..limits.path_conditions import InitialGuess, InitialGuessList
-from ..misc.enums import ControlType, CostType, Shooting
+from ..misc.enums import ControlType, CostType, Shooting, InterpolationType
 from ..misc.utils import check_version
 from ..limits.phase_transition import PhaseTransitionFunctions
 from ..optimization.non_linear_program import NonLinearProgram
@@ -262,7 +262,8 @@ def init_from_initial_guess(sol: list):
             self.vector = np.ndarray((0, 1))
             sol_states, sol_controls = sol[0], sol[1]
             for p, s in enumerate(sol_states):
-                s.init.check_and_adjust_dimensions(self.ocp.nlp[p].nx, self.ocp.nlp[p].ns + 1, "states")
+                ns = self.ocp.nlp[p].ns + 1 if s.init.type != InterpolationType.EACH_FRAME else self.ocp.nlp[p].ns
+                s.init.check_and_adjust_dimensions(self.ocp.nlp[p].nx, ns, "states")
                 for i in range(self.ns[p] + 1):
                     self.vector = np.concatenate((self.vector, s.init.evaluate_at(i)[:, np.newaxis]))
             for p, s in enumerate(sol_controls):
diff --git a/examples/getting_started/example_simulation.py b/examples/getting_started/example_simulation.py
index 1f335de54..6602f8c64 100644
--- a/examples/getting_started/example_simulation.py
+++ b/examples/getting_started/example_simulation.py
@@ -9,8 +9,8 @@
 can be held in the solution. Another goal would be to reload fast a previously saved optimized solution
 """
 
-from bioptim import InitialGuess, Solution, Shooting
-
+from bioptim import InitialGuess, Solution, Shooting, InterpolationType
+import numpy as np
 import pendulum
 
 
@@ -23,6 +23,7 @@
     )
 
     # Simulation the Initial Guess
+    # Interpolation: Constant
     X = InitialGuess([0, 0, 0, 0])
     U = InitialGuess([-1, 1])
 
@@ -32,17 +33,28 @@
     # Uncomment the next line to animate the integration
     # s.animate()
 
+    # Interpolation: Each frame (for instance, values from a previous optimization or from measured data)
+    U = np.random.rand(2, 11)
+    U = InitialGuess(U, interpolation=InterpolationType.EACH_FRAME)
+
+    sol_from_initial_guess = Solution(ocp, [X, U])
+    s = sol_from_initial_guess.integrate(shooting_type=Shooting.SINGLE_CONTINUOUS)
+    print(f"Final position of q from single shooting of initial guess = {s.states['q'][:, -1]}")
+    # Uncomment the next line to animate the integration
+    # s.animate()
+
     # Uncomment the following lines to graph the solution from initial guesses
     # sol_from_initial_guess.graphs(shooting_type=Shooting.SINGLE_CONTINUOUS)
     # sol_from_initial_guess.graphs(shooting_type=Shooting.MULTIPLE)
 
-    # Simulation of the solution. It is not the graph of the solution, it is the graph of a Runge Kutta from the solution
+    # Simulation of the solution. It is not the graph of the solution,
+    # it is the graph of a Runge Kutta from the solution
     sol = ocp.solve()
     s_single = sol.integrate(shooting_type=Shooting.SINGLE_CONTINUOUS)
     # Uncomment the next line to animate the integration
     # s_single.animate()
     print(f"Final position of q from single shooting of the solution = {s_single.states['q'][:, -1]}")
-    s_multiple = sol.integrate(shooting_type=Shooting.MULTIPLE)
+    s_multiple = sol.integrate(shooting_type=Shooting.MULTIPLE, keepdims=False)
     print(f"Final position of q from multiple shooting of the solution = {s_multiple.states['q'][:, -1]}")
 
     # Uncomment the following lines to graph the solution from the actual solution