crtbp_derivs_with_stm Subroutine

    subroutine crtbp_derivs_with_stm(mu,x,dx)

    implicit none

    real(wp),intent(in)                :: mu   !! CRTBP parameter (See [[compute_crtpb_parameter]])
    real(wp),dimension(42),intent(in)  :: x    !! normalized state and STM \([\mathbf{r},\mathbf{v},\mathbf{\Phi}]\)
    real(wp),dimension(42),intent(out) :: dx   !! normalized state and STM derivative \([\dot{\mathbf{r}},\dot{\mathbf{v}},\dot{\mathbf{\Phi}}]\)

    !local variables:
    real,dimension(3) :: rb1,rb2,r,v,g
    real(wp),dimension(6,6) :: A, phi, phi_dot
    real(wp) :: r1,r2,r13,r23,r15,r25,omm,tmu,tomm,c1,c2
    real(wp) :: Uxx,Uxy,Uxz,Uyx,Uyy,Uyz,Uzx,Uzy,Uzz

    !extract variables from x vector:
    r = x(1:3) ! position
    v = x(4:6) ! velocity

    !other parameters:
    omm  = one - mu
    rb1  = [-mu,zero,zero] ! location of body 1
    rb2  = [omm,zero,zero] ! location of body 2
    r1   = norm2(r - rb1)  ! body1 -> sc distance
    r2   = norm2(r - rb2)  ! body2 -> sc distance
    r13  = r1**3
    r23  = r2**3
    r15  = r1**5
    r25  = r2**5
    c1   = omm/r13
    c2   = mu/r23
    tmu  = three*mu
    tomm = three*omm

    !normalized gravity from both bodies:
    g(1) = -c1*(r(1) + mu) - c2*(r(1)-one+mu)
    g(2) = -c1*r(2)        - c2*r(2)
    g(3) = -c1*r(3)        - c2*r(3)

    !STM terms:
    Uxx = one - c1 - c2 + tomm*(r(1)+mu)**2/r15 + tmu*(r(1)-one+mu)**2/r25
    Uyy = one - c1 - c2 + tomm*r(2)**2/r15 + tmu*r(2)**2/r25
    Uzz =       c1 - c2 + tomm*r(3)**2/r15 + tmu*r(3)**2/r25
    Uxy = tomm*(r(1)+mu)*r(2)/r15 + tmu*(r(1)-one+mu)*r(2)/r25
    Uxz = tomm*(r(1)+mu)*r(3)/r15 + tmu*(r(1)-one+mu)*r(3)/r25
    Uyz = tomm* r(2)*r(3)/r15 + tmu*r(2)*r(3)/r25
    Uyx = Uxy
    Uzx = Uxz
    Uzy = Uyz

    !columns of A matrix:
    A(:,1) = [zero,zero,zero,Uxx,Uyx,Uzx]
    A(:,2) = [zero,zero,zero,Uxy,Uyy,Uzy]
    A(:,3) = [zero,zero,zero,Uxz,Uyz,Uzz]
    A(:,4) = [one,zero,zero,zero,-two,zero]
    A(:,5) = [zero,one,zero,two,zero,zero]
    A(:,6) = [zero,zero,one,zero,zero,zero]

    !unpack phi into matrix:
    phi = reshape(x(7:42), shape=[6,6])

    !derivative of phi matrix:
    phi_dot = matmul(A,phi)

    !derivative of x vector:
    dx(1:3)  = v                           ! r_dot
    dx(4)    =  two*v(2) + r(1) + g(1)     ! v_dot
    dx(5)    = -two*v(1) + r(2) + g(2)     !
    dx(6)    =                    g(3)     !
    dx(7:42) = pack (phi_dot, mask=.true.) ! phi_dot

    end subroutine crtbp_derivs_with_stm