anchor.diff:        end do
         do i = 1,NTime
          ATemp = A
  -       RHS = MATMUL(B,ChiTime(:,i)) + MATMUL(D, DPhiW(:,i))
  -       CALL DGESV(N, 1, ATEMP, N, IPIV, RHS, N, INFO)
  -       IF (INFO .NE. 0) THEN
  -         WRITE(*,*) "ERROR IN LINEAR SYSTEM"
  -       END IF
  +       RHS = matmul(B,ChiTime(:,i)) + matmul(D, DPhiW(:,i))
  +       call DGESV(N, 1, ATEMP, N, IPIV, RHS, N, INFO)
  +       if (INFO .NE. 0) then
  +         write(*,*) "ERROR IN LINEAR SYSTEM"
  +       end if
          PhiTime(:,i) = RHS
  -       CALL Wake(N, NWake, NTime, i, PhiTime, DPhiW)
  +       call Wake(N, NWake, NTime, i, PhiTime, DPhiW)
         end do
  -END SUBROUTINE
  +end subroutine
   
   
  -SUBROUTINE SolvePhiLap(N, B, C, DRS, ChiLap, PhiLap)
  +subroutine SolvePhiLap(N, B, C, DRS, ChiLap, PhiLap)
         IMPLICIT NONE
         integer, intent(IN) :: N
         real(kind=8), dimension(N,N), intent(IN) :: B, C
  
positive.code:       do i = 1, N
         A(i,i) = A(i,i) + dcmplx(0.5,0.)
        end do
        RHS = MATMUL(dcmplx(B),ChiLap)
        CALL ZGESV(N, 1, A, N, IPIV, RHS, N, INFO)
        IF (INFO .NE. 0) THEN
          WRITE(*,*) "ERROR IN LINEAR SYSTEM"
        END IF
        PhiLap = RHS
  END SUBROUTINE
  
  
  !Calculate the potential in the field
  !  PHISRF   On surface                          (in)
  !  CHISRF   On surface                          (in)
  !  B, C     In the field                        (in)
  !  PHIFLD   In the field                        (out)
  SUBROUTINE CALCPHIFLD(N, PHISRF, CHISRF, NX, B, C, PHIFLD)
        IMPLICIT NONE
        INTEGER, INTENT(IN) :: N, NX
        REAL(KIND=8), DIMENSION(N), INTENT(IN) :: PHISRF, CHISRF
        REAL(KIND=8), DIMENSION(NX,N) :: B, C
        REAL(KIND=8), DIMENSION(NX), INTENT(OUT) :: PHIFLD
        PHIFLD = MATMUL(B, CHISRF) + MATMUL(C, PHISRF)
  END SUBROUTINE
  
negative.code: subroutine BCondRotLap(Nelem, Xnode, alphaAmpl, ChiLap)
        IMPLICIT NONE
        real(KIND=8), parameter :: PI = 4.D0*datan(1.D0)
        integer, intent(IN) :: Nelem
        real(kind=8), dimension(Nelem,2), intent(IN) :: Xnode
        real(kind=8), intent(IN) :: alphaAmpl
        real(kind=8) :: alpharad
        complex(kind=8), dimension(Nelem), intent(OUT) :: ChiLap
        real(kind=8), dimension(Nelem,2) :: n
        integer :: I
        call normals(Nelem, Xnode, n)
  