anchor.diff:     m%uFlat(10:12) = u%PtfmMesh%RotationVel   (:,1)
      m%uFlat(13:15) = u%PtfmMesh%TranslationAcc(:,1)
      m%uFlat(16:18) = u%PtfmMesh%RotationAcc   (:,1)
  -   !write(*,*)'x%qm',x%qm
  -   !write(*,*)'x%qmdot',x%qmdot
  -   !write(*,*)'m%uFlat',m%uFlat
  -   !write(*,*)'m%F_at_t',m%F_at_t(6+1:6+p%nCB)
  -   !write(*,*)'K22',p%K22
  -   !write(*,*)'C22',p%C22
  -   !write(*,*)'M21',p%M21
  -   !write(*,*)'C21',p%C21
  -
  -
  -   dxdt%qm= x%qmdot
  +
  +   ! --- Computation of qm and qmdot
  +   ! >>> Latex formulae:
      ! \ddot{x2} = -K22 x2 - C22 \dot{x2}  - C21 \dot{x1} - M21 \ddot{x1} + fr2
  -   dxdt%qmdot = - matmul(p%K22,x%qm) - matmul(p%C22,x%qmdot) &
  -                - matmul(p%C21,m%uFlat(7:12)) - matmul(p%M21, m%uFlat(13:18)) + m%F_at_t(6+1:6+p%nCB)
  +   ! >>> MATMUL IMPLEMENTATION 
  +   !dxdt%qm= x%qmdot
  +   !dxdt%qmdot = - matmul(p%K22,x%qm) - matmul(p%C22,x%qmdot) &
  +   !             - matmul(p%C21,m%uFlat(7:12)) - matmul(p%M21, m%uFlat(13:18)) + m%F_at_t(6+1:6+p%nCB)
  +   ! >>> BLAS IMPLEMENTATION 
  +   !           COPY( N   , X                    , INCX, Y      , INCY)
  +   CALL LAPACK_COPY(p%nCB, x%qmdot              , 1  , dxdt%qm    , 1  ) ! qmdot=qmdot
  +   CALL LAPACK_COPY(p%nCB, m%F_at_t(6+1:6+p%nCB), 1  , dxdt%qmdot , 1  )                                          ! qmddot = fr2
  +   !           GEMV(TRS, M    , N     , alpha    , A    , LDA  , X              ,INCX, Beta   ,  Y        , IncY)
  +   CALL LAPACK_GEMV('n', p%nCB, p%nCB , -1.0_ReKi, p%K22, p%nCB, x%qm          , 1  , 1.0_ReKi, dxdt%qmdot, 1   ) !        - K22 x2
  +   CALL LAPACK_GEMV('n', p%nCB, 6     , -1.0_ReKi, p%C21, p%nCB, m%uFlat(7:12) , 1  , 1.0_ReKi, dxdt%qmdot, 1   ) !        - C21 \dot{x1}
  +   CALL LAPACK_GEMV('n', p%nCB, p%nCB , -1.0_ReKi, p%C22, p%nCB, x%qmdot       , 1  , 1.0_ReKi, dxdt%qmdot, 1   ) !        - C22 \dot{x2}
  +   CALL LAPACK_GEMV('n', p%nCB, 6     , -1.0_ReKi, p%M21, p%nCB, m%uFlat(13:18), 1  , 1.0_ReKi, dxdt%qmdot, 1   ) !        - M21 \ddot{x1}
  +
   CONTAINS
       logical function Failed()
           CALL SetErrStatSimple(ErrStat, ErrMsg, 'ExtPtfm_CalcContStateDeriv')
  
positive.code:    m%uFlat(10:12) = u%PtfmMesh%RotationVel   (:,1)
     m%uFlat(13:15) = u%PtfmMesh%TranslationAcc(:,1)
     m%uFlat(16:18) = u%PtfmMesh%RotationAcc   (:,1)
     if (p%nCB>0) then
         ! x flat
         m%xFlat(      1:p%nCB  ) = x%qm   (1:p%nCB)
         m%xFlat(p%nCB+1:2*p%nCB) = x%qmdot(1:p%nCB)
         ! y = Cx + Du + Fy  - TODO consider using lapack for efficiency
         Fc = matmul(p%CMat, m%xFlat) + matmul(p%DMat, m%uFlat) + m%F_at_t(1:6) - matmul(p%M12, m%F_at_t(6+1:6+p%nCB))
     else
         Fc =                           matmul(p%DMat, m%uFlat) + m%F_at_t(1:6) 
     endif
  
negative.code:    TYPE(ExtPtfm_InputType),             INTENT(INOUT) :: u(:)        !< Inputs at t
     REAL(DbKi),                          INTENT(IN   ) :: utimes(:)   !< times of input
     TYPE(ExtPtfm_ParameterType),         INTENT(IN   ) :: p           !< Parameters
     TYPE(ExtPtfm_ContinuousStateType),   INTENT(INOUT) :: x           !< Continuous states at t on input at t + dt on output ! TODO TODO TODO IN
     TYPE(ExtPtfm_DiscreteStateType),     INTENT(IN   ) :: xd          !< Discrete states at t
     TYPE(ExtPtfm_ConstraintStateType),   INTENT(IN   ) :: z           !< Constraint states at t (possibly a guess)
     TYPE(ExtPtfm_OtherStateType),        INTENT(INOUT) :: OtherState  !< Other states at t on input at t + dt on output
  