anchor.diff:  
     write(*,*) "Setting up the test problem"
   
  +  ! ******************************************************
  +  ! (1) create a pendulum body
  +  ! ******************************************************
  +
     call create_body()
  -
  -  print*, alpha
  -  call jac
  -  stop
  +  
  +  ! ******************************************************
  +  ! (2) set the res and jacobian for the linear system
  +  ! ******************************************************
   
     ! set number of equations to solve
     NEQ=12
   
  -  IDID=0  !?
  +  IDID=0  
   
     ! set tolerances
  -  Rtol = 0.005
  -  Atol = 0.005
  +  Rtol = 0.05
  +  Atol = 0.05
   
     ! intial time
     t=0.00
   
     ! final time
     tout=0.03
  -
  -
  -  ! initial state
   
     ! info block
     do ii=1,15
  
positive.code:   info(2)=1
    info(3)=1
  
    ! allocate size of work arrays
    liw=1000
    lrw=1055 
  
  
    ! (1) create a pendulum body
  
    ! (2) set the res and jacobian for the linear system
  
    ! (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)
    end do
  
  end program pendulum
  
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
  