diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py
index a7e0374b2..f4659a423 100644
--- a/rocketpy/simulation/flight.py
+++ b/rocketpy/simulation/flight.py
@@ -502,6 +502,7 @@ def __init__(  # pylint: disable=too-many-arguments,too-many-statements
         verbose=False,
         name="Flight",
         equations_of_motion="standard",
+        use_udot_rail2=True,
         ode_solver="LSODA",
         simulation_mode="6 DOF",
         weathercock_coeff=0.0,
@@ -535,6 +536,11 @@ def __init__(  # pylint: disable=too-many-arguments,too-many-statements
                     x_init, y_init, z_init,
                     vx_init, vy_init, vz_init,
                     e0_init, e1_init, e2_init, e3_init,
+        use_udot_rail2 : bool, optional
+            If True, enable the intermediate 3-DOF rail phase `udot_rail2`.
+            If False, the flight remains in the rail1 phase until the
+            upper rail button exit and then transitions directly to the
+            generalized 6-DOF dynamics. Default is True.
                     w1_init, w2_init, w3_init
                 ]
 
@@ -627,6 +633,8 @@ def __init__(  # pylint: disable=too-many-arguments,too-many-statements
         self.equations_of_motion = equations_of_motion
         self.simulation_mode = simulation_mode
         self.ode_solver = ode_solver
+        # Enable or disable the intermediate 3-DOF rail phase
+        self.use_udot_rail2 = use_udot_rail2
         self.weathercock_coeff = weathercock_coeff
 
         # Controller initialization
@@ -642,6 +650,7 @@ def __init__(  # pylint: disable=too-many-arguments,too-many-statements
         self.flight_phases.add_phase(
             self.t_initial, self.initial_derivative, clear=False
         )
+        # self.flight_phases.add_phase()
         self.flight_phases.add_phase(self.max_time)
 
         # Simulate flight
@@ -785,6 +794,10 @@ def __simulate(self, verbose):
                     break  # Stop simulation if parachute is deployed
 
                 # Step through simulation
+                # Use a local insertion offset `i` for any phases added while
+                # stepping the solver so multiple inserts in the same loop
+                # are placed deterministically (parachute-style pattern).
+                i = 1
                 while phase.solver.status == "running":
                     # Execute solver step, log solution and function evaluations
                     phase.solver.step()
@@ -797,6 +810,192 @@ def __simulate(self, verbose):
                     if verbose:
                         print(f"Current Simulation Time: {self.t:3.4f} s", end="\r")
 
+                    # check for the first time the rocket is between the two rail buttons for tip-off analysis
+                    if len(self.between_rails_state) == 1 and (
+                        self.y_sol[0] ** 2
+                        + self.y_sol[1] ** 2
+                        + (self.y_sol[2] - self.env.elevation) ** 2
+                        >= self.effective_2rl**2
+                    ):
+                        self.between_rails_time = self.t
+                        self.between_rails_time_index = len(self.solution) - 1
+                        self.between_rails_state = self.y_sol
+                        # Optionally insert the udot_rail2 3-DOF phase. If
+                        # disabled, the solver remains in the current phase
+                        # and will transition to generalized dynamics at the
+                        # upper button exit as before.
+                        if self.use_udot_rail2:
+                            self.flight_phases.add_phase(
+                                self.t,
+                                self.udot_rail2,
+                                index=phase_index + i,
+                            )
+                            i += 1
+                            # Prepare to leave loops and start new flight phase
+                            phase.time_nodes.flush_after(node_index)
+                            phase.time_nodes.add_node(self.t, [], [], [])
+                            phase.solver.status = "finished"
+                    # Check for first out of rail event
+                    elif len(self.out_of_rail_state) == 1 and (
+                        self.y_sol[0] ** 2
+                        + self.y_sol[1] ** 2
+                        + (self.y_sol[2] - self.env.elevation) ** 2
+                        >= self.effective_1rl**2
+                    ):
+                        # Check exactly when it went out using root finding
+                        # Disconsider elevation
+                        self.solution[-2][3] -= self.env.elevation
+                        self.solution[-1][3] -= self.env.elevation
+                        # Get points
+                        y0 = (
+                            sum(self.solution[-2][i] ** 2 for i in [1, 2, 3])
+                            - self.effective_1rl**2
+                        )
+                        yp0 = 2 * sum(
+                            self.solution[-2][i] * self.solution[-2][i + 3]
+                            for i in [1, 2, 3]
+                        )
+                        t1 = self.solution[-1][0] - self.solution[-2][0]
+                        y1 = (
+                            sum(self.solution[-1][i] ** 2 for i in [1, 2, 3])
+                            - self.effective_1rl**2
+                        )
+                        yp1 = 2 * sum(
+                            self.solution[-1][i] * self.solution[-1][i + 3]
+                            for i in [1, 2, 3]
+                        )
+                        # Put elevation back
+                        self.solution[-2][3] += self.env.elevation
+                        self.solution[-1][3] += self.env.elevation
+                        # Cubic Hermite interpolation (ax**3 + bx**2 + cx + d)
+                        a, b, c, d = calculate_cubic_hermite_coefficients(
+                            0,
+                            float(phase.solver.step_size),
+                            y0,
+                            yp0,
+                            y1,
+                            yp1,
+                        )
+                        a += 1e-5  # TODO: why??
+                        # Find roots
+                        t_roots = find_roots_cubic_function(a, b, c, d)
+                        # Find correct root
+                        valid_t_root = [
+                            t_root.real
+                            for t_root in t_roots
+                            if 0 < t_root.real < t1 and abs(t_root.imag) < 0.001
+                        ]
+                        if len(valid_t_root) > 1:  # pragma: no cover
+                            raise ValueError(
+                                "Multiple roots found when solving for rail exit time."
+                            )
+                        if len(valid_t_root) == 0:  # pragma: no cover
+                            # Numerical issues can cause no root to be accepted
+                            # by the strict filter. Fall back to a safe
+                            # approximation (midpoint of the interval) so the
+                            # simulation can continue; warn the user so this
+                            # event can be investigated if needed.
+                            warnings.warn(
+                                "No valid roots found when solving for rail exit time; "
+                                "using midpoint fallback.",
+                                RuntimeWarning,
+                            )
+                            fallback = max(min(t1 * 0.5, t1 - 1e-12), 1e-12)
+                            valid_t_root = [fallback]
+                        # Determine final state when upper button is going out of rail
+                        self.t = valid_t_root[0] + self.solution[-2][0]
+                        interpolator = phase.solver.dense_output()
+                        self.y_sol = interpolator(self.t)
+                        self.solution[-1] = [self.t, *self.y_sol]
+                        self.out_of_rail_time = self.t
+                        self.out_of_rail_time_index = len(self.solution) - 1
+                        self.out_of_rail_state = self.y_sol
+                        # Create new flight phase
+                        self.flight_phases.add_phase(
+                            self.t,
+                            self.u_dot_generalized,
+                            index=phase_index + i,
+                        )
+                        i += 1
+                        # Prepare to leave loops and start new flight phase
+                        phase.time_nodes.flush_after(node_index)
+                        phase.time_nodes.add_node(self.t, [], [], [])
+                        phase.solver.status = "finished"
+
+                    # Check for apogee event
+                    # TODO: negative vz doesn't really mean apogee. Improve this.
+                    if len(self.apogee_state) == 1 and self.y_sol[5] < 0:
+                        # Assume linear vz(t) to detect when vz = 0
+                        t0, vz0 = self.solution[-2][0], self.solution[-2][6]
+                        t1, vz1 = self.solution[-1][0], self.solution[-1][6]
+                        t_root = find_root_linear_interpolation(t0, t1, vz0, vz1, 0)
+                        # Fetch state at t_root
+                        interpolator = phase.solver.dense_output()
+                        self.apogee_state = interpolator(t_root)
+                        # Store apogee data
+                        self.apogee_time = t_root
+                        self.apogee_x = self.apogee_state[0]
+                        self.apogee_y = self.apogee_state[1]
+                        self.apogee = self.apogee_state[2]
+
+                        if self.terminate_on_apogee:
+                            self.t = self.t_final = t_root
+                            # Roll back solution
+                            self.solution[-1] = [self.t, *self.apogee_state]
+                            # Set last flight phase
+                            self.flight_phases.flush_after(phase_index)
+                            self.flight_phases.add_phase(self.t)
+                            # Prepare to leave loops and start new flight phase
+                            phase.time_nodes.flush_after(node_index)
+                            phase.time_nodes.add_node(self.t, [], [], [])
+                            phase.solver.status = "finished"
+                        elif len(self.solution) > 2:
+                            # adding the apogee state to solution increases accuracy
+                            # we can only do this if the apogee is not the first state
+                            self.solution.insert(-1, [t_root, *self.apogee_state])
+                    # Check for impact event
+                    if self.y_sol[2] < self.env.elevation:
+                        # Check exactly when it happened using root finding
+                        # Cubic Hermite interpolation (ax**3 + bx**2 + cx + d)
+                        a, b, c, d = calculate_cubic_hermite_coefficients(
+                            x0=0,  # t0
+                            x1=float(phase.solver.step_size),  # t1 - t0
+                            y0=float(self.solution[-2][3] - self.env.elevation),  # z0
+                            yp0=float(self.solution[-2][6]),  # vz0
+                            y1=float(self.solution[-1][3] - self.env.elevation),  # z1
+                            yp1=float(self.solution[-1][6]),  # vz1
+                        )
+                        # Find roots
+                        t_roots = find_roots_cubic_function(a, b, c, d)
+                        # Find correct root
+                        t1 = self.solution[-1][0] - self.solution[-2][0]
+                        valid_t_root = [
+                            t_root.real
+                            for t_root in t_roots
+                            if abs(t_root.imag) < 0.001 and 0 < t_root.real < t1
+                        ]
+                        if len(valid_t_root) > 1:  # pragma: no cover
+                            raise ValueError(
+                                "Multiple roots found when solving for impact time."
+                            )
+                        # Determine impact state at t_root
+                        self.t = self.t_final = valid_t_root[0] + self.solution[-2][0]
+                        interpolator = phase.solver.dense_output()
+                        self.y_sol = self.impact_state = interpolator(self.t)
+                        # Roll back solution
+                        self.solution[-1] = [self.t, *self.y_sol]
+                        # Save impact state
+                        self.x_impact = self.impact_state[0]
+                        self.y_impact = self.impact_state[1]
+                        self.z_impact = self.impact_state[2]
+                        self.impact_velocity = self.impact_state[5]
+                        # Set last flight phase
+                        self.flight_phases.flush_after(phase_index)
+                        self.flight_phases.add_phase(self.t)
+                        # Prepare to leave loops and start new flight phase
+                        phase.time_nodes.flush_after(node_index)
+                        phase.time_nodes.add_node(self.t, [], [], [])
+                        phase.solver.status = "finished"
                     if self.__check_simulation_events(phase, phase_index, node_index):
                         break  # Stop if simulation termination event occurred
 
@@ -1459,6 +1658,9 @@ def __init_solution_monitors(self):
         self.out_of_rail_time = 0
         self.out_of_rail_time_index = 0
         self.out_of_rail_state = np.array([0])
+        self.between_rails_state = np.array([0])
+        self.between_rails_time = 0
+        self.between_rails_time_index = 0
         self.apogee_state = np.array([0])
         self.apogee = 0
         self.apogee_time = 0
@@ -1500,6 +1702,18 @@ def __init_flight_state(self):
             e0_init, e1_init, e2_init, e3_init = euler313_to_quaternions(
                 self.phi_init, self.theta_init, self.psi_init
             )
+
+            K_init = Matrix.transformation([e0_init, e1_init, e2_init, e3_init])
+
+            # Body axis pointing along rocket symmetry (body z-axis)
+            body_axis = Vector([0, 0, 1])
+
+            # Attitude vector in inertial frame
+            attitude_vec = K_init @ body_axis
+
+            # Unit vector (normalize)
+            self.attitude_unit = attitude_vec / abs(attitude_vec)
+
             # Store initial conditions
             self.initial_solution = [
                 self.t_initial,
@@ -1527,6 +1741,10 @@ def __init_flight_state(self):
             self.out_of_rail_state = self.initial_solution[1:]
             self.out_of_rail_time = self.initial_solution[0]
             self.out_of_rail_time_index = 0
+            # save out of rail 2 state and time with the same data as out of rail
+            self.between_rails_state = self.initial_solution[1:]
+            self.between_rails_time = self.initial_solution[0]
+            self.between_rails_time_index = 0
             # Set initial derivative for 6-DOF flight phase
             self.initial_derivative = self.u_dot_generalized
         else:
@@ -1535,6 +1753,10 @@ def __init_flight_state(self):
             self.out_of_rail_state = self.initial_solution[1:]
             self.out_of_rail_time = self.initial_solution[0]
             self.out_of_rail_time_index = 0
+            # save out of rail 2 state and time with the same data as out of rail
+            self.between_rails_state = self.initial_solution[1:]
+            self.between_rails_time = self.initial_solution[0]
+            self.between_rails_time_index = 0
             self.t_initial = self.initial_solution[0]
             self.initial_derivative = self.u_dot_generalized
         if self._controllers or self.sensors:
@@ -1874,7 +2096,7 @@ def udot_rail1(self, t, u, post_processing=False):
         return [vx, vy, vz, ax, ay, az, 0, 0, 0, 0, 0, 0, 0]
 
     def udot_rail2(self, t, u, post_processing=False):  # pragma: no cover
-        """[Still not implemented] Calculates derivative of u state vector with
+        """[WIP] Calculates derivative of u state vector with
         respect to time when rocket is flying in 3 DOF motion in the rail.
 
         Parameters
@@ -1894,8 +2116,184 @@ def udot_rail2(self, t, u, post_processing=False):  # pragma: no cover
             State vector defined by u_dot = [vx, vy, vz, ax, ay, az,
             e0dot, e1dot, e2dot, e3dot, alpha1, alpha2, alpha3].
         """
-        # Hey! We will finish this function later, now we just can use u_dot
-        return self.u_dot_generalized(t, u, post_processing=post_processing)
+
+        # Retrieve integration data
+        _, _, z, vx, vy, vz, e0, e1, e2, e3, omega1, omega2, omega3 = u
+
+        # Create necessary vectors
+        # r = Vector([x, y, z])  # CDM position vector
+        v = Vector([vx, vy, vz])  # CDM velocity vector
+        e = [e0, e1, e2, e3]  # Euler parameters/quaternions
+        w = Vector([omega1, omega2, omega3])  # Angular velocity vector
+
+        # Retrieve necessary quantities
+        ## Rocket mass
+        total_mass = self.rocket.total_mass.get_value_opt(t)
+        total_mass_dot = self.rocket.total_mass_flow_rate.get_value_opt(t)
+        total_mass_ddot = self.rocket.total_mass_flow_rate.differentiate_complex_step(t)
+        ## CM position vector and time derivatives relative to CDM in body frame
+        r_CM_z = self.rocket.com_to_cdm_function
+        r_CM_t = r_CM_z.get_value_opt(t)
+        r_CM = Vector([0, 0, r_CM_t])
+        r_CM_dot = Vector([0, 0, r_CM_z.differentiate_complex_step(t)])
+        r_CM_ddot = Vector([0, 0, r_CM_z.differentiate(t, order=2)])
+        ## Nozzle position vector
+        r_NOZ = Vector([0, 0, self.rocket.nozzle_to_cdm])
+        ## Nozzle gyration tensor
+        S_nozzle = self.rocket.nozzle_gyration_tensor
+        ## Inertia tensor
+        inertia_tensor = self.rocket.get_inertia_tensor_at_time(t)
+        ## Inertia tensor time derivative in the body frame
+        I_dot = self.rocket.get_inertia_tensor_derivative_at_time(t)
+
+        # Calculate the Inertia tensor relative to CM
+        H = (r_CM.cross_matrix @ -r_CM.cross_matrix) * total_mass
+        I_CM = inertia_tensor - H
+
+        # Prepare transformation matrices
+        K = Matrix.transformation(e)
+        Kt = K.transpose
+
+        # Compute aerodynamic forces and moments
+        R1, R2, R3, M1, M2, M3 = 0, 0, 0, 0, 0, 0
+
+        ## Drag force
+        rho = self.env.density.get_value_opt(z)
+        wind_velocity_x = self.env.wind_velocity_x.get_value_opt(z)
+        wind_velocity_y = self.env.wind_velocity_y.get_value_opt(z)
+        wind_velocity = Vector([wind_velocity_x, wind_velocity_y, 0])
+        free_stream_speed = abs((wind_velocity - Vector(v)))
+        speed_of_sound = self.env.speed_of_sound.get_value_opt(z)
+        free_stream_mach = free_stream_speed / speed_of_sound
+
+        if self.rocket.motor.burn_start_time < t < self.rocket.motor.burn_out_time:
+            pressure = self.env.pressure.get_value_opt(z)
+            net_thrust = max(
+                self.rocket.motor.thrust.get_value_opt(t)
+                + self.rocket.motor.pressure_thrust(pressure),
+                0,
+            )
+            drag_coeff = self.rocket.power_on_drag.get_value_opt(free_stream_mach)
+        else:
+            net_thrust = 0
+            drag_coeff = self.rocket.power_off_drag.get_value_opt(free_stream_mach)
+        R3 += -0.5 * rho * (free_stream_speed**2) * self.rocket.area * drag_coeff
+        # Get rocket velocity in body frame
+        velocity_in_body_frame = Kt @ v
+        # Calculate lift and moment for each component of the rocket
+        for aero_surface, _ in self.rocket.aerodynamic_surfaces:
+            # Component cp relative to CDM in body frame
+            comp_cp = self.rocket.surfaces_cp_to_cdm[aero_surface]
+            # Component absolute velocity in body frame
+            comp_vb = velocity_in_body_frame + (w ^ comp_cp)
+            # Wind velocity at component altitude
+            comp_z = z + (K @ comp_cp).z
+            comp_wind_vx = self.env.wind_velocity_x.get_value_opt(comp_z)
+            comp_wind_vy = self.env.wind_velocity_y.get_value_opt(comp_z)
+            # Component freestream velocity in body frame
+            comp_wind_vb = Kt @ Vector([comp_wind_vx, comp_wind_vy, 0])
+            comp_stream_velocity = comp_wind_vb - comp_vb
+            comp_stream_speed = abs(comp_stream_velocity)
+            comp_stream_mach = comp_stream_speed / speed_of_sound
+            # Reynolds at component altitude
+            # TODO: Reynolds is only used in generic surfaces. This calculation
+            # should be moved to the surface class for efficiency
+            comp_reynolds = (
+                self.env.density.get_value_opt(comp_z)
+                * comp_stream_speed
+                * aero_surface.reference_length
+                / self.env.dynamic_viscosity.get_value_opt(comp_z)
+            )
+            # Forces and moments
+            X, Y, Z, M, N, L = aero_surface.compute_forces_and_moments(
+                comp_stream_velocity,
+                comp_stream_speed,
+                comp_stream_mach,
+                rho,
+                comp_cp,
+                w,
+                comp_reynolds,
+            )
+            R1 += X
+            R2 += Y
+            R3 += Z
+            M1 += M
+            M2 += N
+            M3 += L
+
+        # Off center moment
+        M1 += (
+            self.rocket.cp_eccentricity_y * R3
+            + self.rocket.thrust_eccentricity_y * net_thrust
+        )
+        M2 -= (
+            self.rocket.cp_eccentricity_x * R3
+            + self.rocket.thrust_eccentricity_x * net_thrust
+        )
+        M3 += self.rocket.cp_eccentricity_x * R2 - self.rocket.cp_eccentricity_y * R1
+
+        weight_in_body_frame = Kt @ Vector(
+            [0, 0, -total_mass * self.env.gravity.get_value_opt(z)]
+        )
+
+        T00 = total_mass * r_CM
+        T03 = 2 * total_mass_dot * (r_NOZ - r_CM) - 2 * total_mass * r_CM_dot
+        T04 = (
+            Vector([0, 0, net_thrust])
+            - total_mass * r_CM_ddot
+            - 2 * total_mass_dot * r_CM_dot
+            + total_mass_ddot * (r_NOZ - r_CM)
+        )
+        T05 = total_mass_dot * S_nozzle - I_dot
+
+        T20 = (
+            ((w ^ T00) ^ w)
+            + (w ^ T03)
+            + T04
+            + weight_in_body_frame
+            + Vector([R1, R2, R3])
+        )
+
+        T21 = (
+            ((inertia_tensor @ w) ^ w)
+            + T05 @ w
+            - (weight_in_body_frame ^ r_CM)
+            + Vector([M1, M2, M3])
+        )
+
+        # Angular velocity derivative
+        w_dot = I_CM.inverse @ (T21 + (T20 ^ r_CM))
+        # Enforce zero roll acceleration for 3-DOF rail motion by creating
+        # a new Vector with the third component set to zero instead of
+        # attempting item assignment on the Vector type.
+        w_dot = Vector([w_dot[0], w_dot[1], 0.0])
+
+        # Euler parameters derivative
+        e_dot = [
+            0.5 * (-omega1 * e1 - omega2 * e2),  # - omega3 * e3),
+            0.5 * (omega1 * e0 - omega2 * e3),  # omega3 * e2
+            0.5 * (omega2 * e0 + omega1 * e3),  # - omega3 * e1
+            0.5 * (omega2 * e1 - omega1 * e2),  # +omega3 * e0
+        ]
+
+        # Velocity vector derivative + Coriolis acceleration
+        w_earth = Vector(self.env.earth_rotation_vector)
+        v_dot = K @ (T20 / total_mass - (r_CM ^ w_dot)) - 2 * (w_earth ^ v)
+
+        # Position vector derivative: projection of velocity along the rail
+        rail = self.attitude_unit  # unit vector inertial frame
+        velocity_vec = Vector([vx, vy, vz])
+        r_dot = rail * (velocity_vec @ rail)
+
+        # Create u_dot
+        u_dot = [*r_dot, *v_dot, *e_dot, *w_dot]
+
+        if post_processing:
+            self.__post_processed_variables.append(
+                [t, *v_dot, *w_dot, R1, R2, R3, M1, M2, M3, net_thrust]
+            )
+
+        return u_dot
 
     def u_dot(self, t, u, post_processing=False):  # pylint: disable=too-many-locals,too-many-statements
         """Calculates derivative of u state vector with respect to time
@@ -4106,6 +4504,7 @@ def to_dict(self, **kwargs):
                     "time": self.time,
                     "out_of_rail_velocity": self.out_of_rail_velocity,
                     "out_of_rail_state": self.out_of_rail_state,
+                    "between_rails_state": self.between_rails_state,
                     "apogee_x": self.apogee_x,
                     "apogee_y": self.apogee_y,
                     "apogee_state": self.apogee_state,
diff --git a/tests/unit/simulation/test_udot_rail2_feature.py b/tests/unit/simulation/test_udot_rail2_feature.py
new file mode 100644
index 000000000..b50122c8c
--- /dev/null
+++ b/tests/unit/simulation/test_udot_rail2_feature.py
@@ -0,0 +1,212 @@
+"""Unit tests for the udot_rail2 rail-phase feature.
+
+These tests follow the project's testing conventions: each test is named
+`test_methodname`, uses the Arrange / Act / Assert pattern, and the expected
+behaviour is documented in the test docstring.
+
+Coverage:
+- Phase insertion ordering when `use_udot_rail2` is enabled
+- udot_rail2 enforces zero roll acceleration and projects `r_dot` on the rail
+- CSV comparison generation for enabled/disabled runs
+
+These tests are intentionally written to avoid plotting and optional-dependency
+features so they run reliably in CI and local environments.
+"""
+
+import csv
+import math
+import os
+
+import numpy as np
+
+from rocketpy.mathutils import Matrix, Vector
+
+
+def _yaw_deg(v: Vector):
+    return math.degrees(math.atan2(v.y, v.x))
+
+
+def _pitch_deg(v: Vector):
+    return math.degrees(math.atan2(v.z, math.hypot(v.x, v.y)))
+
+
+def _body_axis_from_e(e):
+    K = Matrix.transformation(e)
+    return K @ Vector([0.0, 0.0, 1.0])
+
+
+def test_udot_rail2_inserts_phase_in_order(calisto_robust, example_spaceport_env):
+    """When `use_udot_rail2=True`, the intermediate 3-DOF `udot_rail2` phase
+    is inserted before the 6-DOF generalized phase (`u_dot_generalized`).
+
+    Arrange: build a Flight with `use_udot_rail2=True`.
+    Act: inspect `flight.flight_phases` derivatives names.
+    Assert: `udot_rail2` appears before `u_dot_generalized` in the phase list.
+    """
+    # Arrange
+    from rocketpy.simulation.flight import Flight
+
+    flight = Flight(
+        rocket=calisto_robust,
+        environment=example_spaceport_env,
+        rail_length=5.2,
+        inclination=85,
+        heading=0,
+        terminate_on_apogee=False,
+        use_udot_rail2=True,
+    )
+
+    # Act
+    derivative_names = [
+        phase.derivative.__name__ if phase.derivative is not None else None
+        for phase in flight.flight_phases.list
+    ]
+
+    # Assert
+    assert "udot_rail2" in derivative_names, "udot_rail2 phase not present"
+    assert "u_dot_generalized" in derivative_names, (
+        "u_dot_generalized phase not present"
+    )
+    assert derivative_names.index("udot_rail2") < derivative_names.index(
+        "u_dot_generalized"
+    ), "udot_rail2 should be inserted before u_dot_generalized"
+
+
+def test_udot_rail2_no_roll_and_alignment(calisto_robust, example_spaceport_env):
+    """udot_rail2 must enforce zero roll acceleration and set `r_dot` as the
+    projection of the velocity vector onto the inertial rail axis.
+
+    Arrange: create flight with `use_udot_rail2=True` and find the between-rails
+    time/state. Act: evaluate `udot_rail2(t, u)` at that instant. Assert: the
+    angular-acceleration third component (roll) is zero and `r_dot` equals the
+    projection of velocity on `flight.attitude_unit`.
+    """
+    from rocketpy.simulation.flight import Flight
+
+    # Arrange
+    flight = Flight(
+        rocket=calisto_robust,
+        environment=example_spaceport_env,
+        rail_length=5.2,
+        inclination=85,
+        heading=0,
+        terminate_on_apogee=False,
+        use_udot_rail2=True,
+    )
+
+    # Act
+    t_between = getattr(flight, "between_rails_time", None)
+    u_between = getattr(flight, "between_rails_state", None)
+
+    # If the flight never registered between-rails, skip the detailed asserts
+    # (the test still passes as it did not exercise the condition).
+    if t_between is None or u_between is None:
+        return
+
+    u_dot = flight.udot_rail2(t_between, u_between)
+
+    # u_dot layout for udot_rail2: [r_dot_x, r_dot_y, r_dot_z, v_dot_x, v_dot_y, v_dot_z, e_dot..., w_dot_x, w_dot_y, w_dot_z]
+    r_dot = Vector(u_dot[0:3])
+    v_dot = Vector(u_dot[3:6])
+    # angular accelerations are last three entries
+    w_dot = Vector(u_dot[-3:])
+
+    # Assert: roll acceleration is zero (third component)
+    assert abs(w_dot[2]) < 1e-12, f"Expected zero roll acceleration, got {w_dot[2]}"
+
+    # Assert: r_dot is projection of velocity onto rail axis
+    rail = Vector(flight.attitude_unit)
+    velocity = Vector(u_between[3:6])
+    projected = rail * (velocity @ rail)
+
+    diff = r_dot - projected
+    assert float(abs(diff)) < 1e-8, (
+        f"r_dot not a projection onto rail (err={float(abs(diff))})"
+    )
+
+
+def test_udot_rail2_csv_comparison_generation(
+    calisto_robust, example_spaceport_env, tmp_path
+):
+    """Generate CSVs comparing runs with udot_rail2 enabled/disabled.
+
+    Arrange: create two flights (enabled/disabled). Act: write CSVs with
+    between-rails and out-of-rail pitch/yaw. Assert: files exist and contain
+    the expected header row; also return numeric deltas for quick inspection.
+    """
+    from rocketpy.simulation.flight import Flight
+
+    out_dir = tmp_path / "udot_rail2_output"
+    out_dir.mkdir(parents=True, exist_ok=True)
+
+    results = []
+
+    for enabled in (True, False):
+        # Arrange / Act
+        flight = Flight(
+            rocket=calisto_robust,
+            environment=example_spaceport_env,
+            rail_length=5.2,
+            inclination=85,
+            heading=0,
+            terminate_on_apogee=False,
+            use_udot_rail2=enabled,
+        )
+
+        t_between = getattr(flight, "between_rails_time", None)
+        t_out = getattr(flight, "out_of_rail_time", None)
+
+        def sample_at(t):
+            if t is None:
+                return None, None
+            sol = min(flight.solution, key=lambda row: abs(row[0] - t))
+            e = sol[7:11]
+            body = _body_axis_from_e(e)
+            return _pitch_deg(body), _yaw_deg(body)
+
+        pitch_between, yaw_between = sample_at(t_between)
+        pitch_out, yaw_out = sample_at(t_out)
+
+        tag = "enabled" if enabled else "disabled"
+        csv_path = out_dir / f"calisto_angles_udot_rail2_{tag}.csv"
+        with open(csv_path, "w", newline="") as f:
+            writer = csv.writer(f)
+            writer.writerow(["time_event", "pitch_deg", "yaw_deg"])
+            if t_between is not None:
+                writer.writerow(["between_rails", pitch_between, yaw_between])
+            if t_out is not None:
+                writer.writerow(["out_of_rail", pitch_out, yaw_out])
+
+        # Assert: file created and header present
+        assert csv_path.exists(), f"CSV was not created: {csv_path}"
+        with open(csv_path, "r", newline="") as f:
+            lines = f.read().splitlines()
+        assert lines and lines[0].startswith("time_event,pitch_deg,yaw_deg"), (
+            "CSV header mismatch"
+        )
+
+        results.append(
+            (
+                enabled,
+                t_between,
+                pitch_between,
+                yaw_between,
+                t_out,
+                pitch_out,
+                yaw_out,
+                str(csv_path),
+            )
+        )
+
+    # Provide a final sanity check: both CSVs were created
+    assert all(os.path.exists(r[-1]) for r in results)
+
+    # Optional: compute numeric deltas for out_of_rail if both present
+    enabled_row = next(r for r in results if r[0] is True)
+    disabled_row = next(r for r in results if r[0] is False)
+
+    if enabled_row[4] is not None and disabled_row[4] is not None:
+        delta_pitch = abs((enabled_row[5] or 0) - (disabled_row[5] or 0))
+        delta_yaw = abs((enabled_row[6] or 0) - (disabled_row[6] or 0))
+        # The test only asserts that deltas are numeric and finite
+        assert math.isfinite(delta_pitch) and math.isfinite(delta_yaw)