anchor.diff:          ! (c) Iterative solution (should implement)
           !-------------------------------------------------------------!
   
  -        !jac = direct_solve(jac); dq = matmul(jac, -res);
  -        
  -        dq = direct_solve(jac, -res, TOT_NDOF)
  -        
  +        jac = direct_solve(jac); dq = matmul(jac, -res);
  +
  +        !dq = direct_solve(jac, -res, TOT_NDOF)
  +
           !-------------------------------------------------------------!
           ! Calcualte residual tolratances                    
           !-------------------------------------------------------------!
   
           update_norm  =  maxval(abs(dq)); res_norm  =  maxval(abs(res));
   
  -        if (k .eq. 1)   write(filenum,'(i4, 25F25.16)') newton_cnt, &
  -             &(q(i), i = 1, size(q)),  (q_dot(i),i = 1, size(q_dot)), update_norm
  +        if (k .eq. 1)   write(filenum,'(i4, 26F25.16)') &
  +             & newton_cnt, &
  +             & update_norm, res_norm, &
  +             & (q(i), i = 1, size(q)), &
  +             & (q_dot(i),i = 1, size(q_dot))
   
           !-------------------------------------------------------------!
           ! Update the state variables and then update the body
  
positive.code: end program pendulum
  
  !*********************************************************************!
  ! Callbaclk function for approximation of jacobian using finite-diff
  !*********************************************************************!
  subroutine residual2(q, qdot, f)
  
  
negative.code: 
          q            = get_updated_q(q, dq)
          q_dot        = get_updated_q_dot(q_dot, dq)
          !q_double_dot = get_updated_q_double_dot(q_double_dot, dq) 
  
          !-------------------------------------------------------------!
          ! Check for convergence and stop is necessary
          !-------------------------------------------------------------!
  
          !is converged?
          if ( update_norm .le. ABS_TOL .AND. res_norm .le. ABS_TOL) then
             
             exit newton
          end if
  
          ! is disverged?
          if ( update_norm .ge. 1.0d5 .AND. res_norm .le. 1.0d5) then           
             call disp(" > Solution diverging - aborting time integrtn")
             exit newton
          end if
  
          !? max iters reached
          if (k .eq. MAX_NEWTON_ITER) then
             call disp(" >> Newton solution failed in ", k , "iterations")
  