anchor.diff:       newton: do k = 1, MAX_NEWTON_ITER
   
           ! (3) Assemble the residual and jacobian using the body
  -        
  +
           ! call this method to update the state for every other iter
           if (k .gt. 1) then
              call set_state(q, q_dot, body1)
           end if
   
           ! residual assembly
  -        call disp(" >> Assembling the residuals...")
  +        !call disp(" >> Assembling the residuals...")
           res  = get_residual(body1)
  -        call disp(" >> Residual assembly complete...")
  -        call disp("   R   =   ", res)
  +        !call disp(" >> Residual assembly complete...")
  +        !call disp("   R   =   ", res)
   
           ! jacobian assembly
  -        call disp(" >> Assembling the Jacobian...")
  +        !call disp(" >> Assembling the Jacobian...")
           jac  = get_jacobian(body1)
  -        call disp(" >> Jacobian assembly complete...")
  -        call disp("   JAC =   ",jac)
  +        !call disp(" >> Jacobian assembly complete...")
  +        !call disp("   JAC =   ",jac)
   
           ! ******************************************************
           ! (3) Solve the linear system and compute update delta_q
           ! ******************************************************
   
  -        call disp(" >> Calling the linear solver...") 
  +        !call disp(" >> Calling the linear solver...") 
   
           !dq = linear_solve(jac,-res,'GMRES')
   
  -        call disp(" >> Update to residual is computed...") 
  +        !call disp(" >> Update to residual is computed...") 
   
           l2_norm   = norm2(dq); linf_norm =  maxval(abs(dq));
   
  
positive.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)
  
          call disp("   q      =   ", q)
  !        call disp("   qdot   =   ", q_dot)
  
          if (l2_norm .le. REL_TOL .OR. linf_norm.le. ABS_TOL) &
               & exit newton
  
negative.code: 
    call disp(" >> End of program...") 
    call disp("========================================================")
  
  end program pendulum
  
  !*********************************************************************!
  ! Callbaclk function for approximation of jacobian using finite-diff
  !*********************************************************************!
  subroutine residual2(q, qdot, f)
  
    use global_constants, only: sp, dp, TOT_NDOF, NUM_SPAT_DIM
    use types
    use bodies
    use residual
  
    real(dp), intent(in), dimension(TOT_NDOF)  :: q, qdot
    real(dp), intent(out), dimension(TOT_NDOF) :: f
    real(dp), dimension(NUM_SPAT_DIM)  :: re    
    real(dp)   :: mass
    type(body) :: body1
  
    mass  = 2.0_dp;   re   = (/ 0.1_dp, 0.2_dp, 0.3_dp /);
    body1 = create_body(mass, vector(re), q, qdot);
  
    f = get_residual(body1);  
  
  !!$  f = (/ sin(q(1))*sin(q(3)), sin(q(2)), sin(q(3)), sin(q(4)),&
  !!$       & sin(q(5)), sin(q(6)), &
  !!$       &sin(q(7)), sin(q(8)), sin(q(9)), &
  !!$       &sin(q(10)), sin(q(11)), sin(q(12)*sin(q(1))) /)
  !!$
  
  end subroutine residual2
  
  !*********************************************************************!
  ! Callbaclk function for approximation of jacobian using finite-diff
  !*********************************************************************!
  subroutine residual(q, f)
  
    use global_constants, only: sp, dp, TOT_NDOF
  