anchor.diff:    info(2)=1
     info(3)=1
   
  +  ! derivatives automatically by numerical differences
  +  info(5) = 1   ! Yes - Set INFO(5)=0   ! No  - Set INFO(5)=1 ! provide JAC
  +
  +  info(11) = 0 ! 1 = Y and Y prime are consistent ; 0 = not consistent compute automatically
  +
     ! allocate size of work arrays
  -  liw=1000
  -  lrw=1055 
  -
  -
  -  ! (1) create a pendulum body
  -
  -  ! (2) set the res and jacobian for the linear system
  +  liw=1000;  lrw=1055 ;
   
     ! (3) call the solver
  -  write(*,*) '***************'
  -  write(*,*) 'time,                              positions     ' 
  -  do  while (t .le. tout) 
  -     call ddassl(res,neq,t,y,yprime,tout,info,rtol,atol, idid,rwork,lrw,iwork,liw,rpar,ipar,jac)
  -     write(*,'(8f13.8)') t,y(1),y(2),y(3),y(4),y(5),y(6),y(7)
  +  write(*,*) '************************************************************************************&
  +       &*************************************************************************************'
  +  write(*,*) '       Time', '            q1', '           q2', '          q3', '            q4', &
  +       &'          q5', '            q6', '           q7', '          q8',&
  +       & '          q9', '           q10', '         q11',  '           q12'
  +  write(*,*) '************************************************************************************&
  +       &*************************************************************************************'
  +  do while (t .le. tout) 
  +     call ddassl(res,neq,t,y,yprime,tout,info,rtol,atol,&
  +          & idid,rwork,lrw,iwork,liw,rpar,ipar,jac)
  +     write(*,'(8f13.2)') t,y(1),y(2),y(3),y(4),y(5),y(6),y(7),y(8),y(9),y(10),y(11),y(12)
     end do
   
   end program pendulum
  
positive.code: !!$
  
    ! define positions of body frame from inertial (q terms)
    r       = (/ 2.0_dp, 2._dp, 2.0_dp /)
    theta   = (/ deg2rad(30.0d0), deg2rad(45.0d0), deg2rad(60.0d0) /)
    v     = (/ 0.0d0, 0.0d0, 0.0d0 /)
    omega = (/ 0.0d0, 0.0d0, 1.0d0 /)
  
    alpha%r         = vector(r)
    alpha%theta     = vector(theta)
  
negative.code:   alpha%gr = alpha%J*alpha%omega_dot
  
  end subroutine create_body
  
  
  
  subroutine disp(t,y)
  
    use constants,only:dp,sp
    implicit none
  
    real(dp)    :: t, y(12)
    integer(sp) :: i
  
    write(*,'(13f13.2)') t, (y(i),i=1,12)
  
  
  end subroutine disp
  
  
  subroutine disp_mat(t, A,m,n)
  
    use constants, only:dp,sp
    implicit none
  
    real(dp), intent(in) :: t
    integer(sp)          :: i, j, m, n
    real(dp)             :: A(m,n)
  
  
    do i = 1, n
       write(*,'(13f13.2)') t, (A(j,i),j=1,m)
    end do
  
  
  end subroutine disp_mat
  