anchor.diff:  
     !?? get the residual
     delta = get_vector_elements(R_rigid(alpha),4)
  -  print*, "residual vector = ", delta
  +!  print*, "residual vector = ", delta
  +
  +  call disp(t,delta)
   
   end subroutine res
   
   !********************************
   ! implementation of the jacobian
   !********************************
  -subroutine jac
  +subroutine jac(t,y,yprime,pd,cj,rpar,ipar)
  +!!$C  to define the matrix of partial derivatives
  +!!$C             PD=DG/DY+CJ*DG/DYPRIME
  +!!$C         CJ is a scalar which is input to JAC.
  +!!$C         For the given values of T,Y,YPRIME, the
  +!!$C         subroutine must evaluate the non-zero partial
  +!!$C         derivatives for each equation and each solution
  +!!$C         component, and store these values in the
  +!!$C         matrix PD.  The elements of PD are set to zero
  +!!$C         before each call to JAC so only non-zero elements
  +!!$C         need to be defined.
  +!!$C
  +!!$C         Subroutine JAC must not alter T,Y,(*),YPRIME(*), or CJ.
  +!!$C         You must declare the name JAC in an EXTERNAL statement in
  +!!$C         your program that calls DDASSL.  You must dimension Y,
  +!!$C         YPRIME and PD in JAC.
  +
     use constants,only:dp,sp
     use body_mod,only:alpha
     use jacobian_mod,only:jac_rigid
  -  real(dp) :: a
  +  implicit none
  +  real(dp)     :: cj
  +  real(dp)     :: t
  +  real(dp)     :: y(12)
  +  real(dp)     :: yprime(12)
  +  integer(sp)  :: ires
  +  real(dp)     :: rpar(*)
  +  integer(sp)  :: ipar(*)
  +  real(dp)     :: pd(12,12)
  +  integer(sp)  :: i,j,k
  +
     ! JACOBIAN  = jac_rigid(alpha,a) !?? implement jacobian
  -  a = 0.1d0
  -
  -  print*, "JAC:", jac_rigid (alpha, a)
  -
  -  stop"where is the impl"
  + ! print*,cj 
  +!  stop
  +  cj = 0.1
  +  pd = jac_rigid (alpha, cj)
  +
  +  print*,  shape(pd)
  +  
  +  call disp_mat(t, pd, 12, 12)
  +
  +  stop
  +
  +  !  stop"where is the impl"
   
   end subroutine jac
   
  
positive.code: 
    implicit none
  
  !  type(body) :: alpha
    real(dp)   :: residual(12)
    real(dp)   :: jacob(12,12)
  
    ! solver related variables
    integer neq,info(15),idid,lrw,liw,iwork(1000),ipar, my, ii
    double precision t, y(12), yprime(12), tout, rtol(12),atol(12), rwork(1055), rpar, h0
  
negative.code:   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 ;
  
    ! (3) call the solver
    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
  