subroutine dglsbp(nu, order, bc, e)
      integer nu
      integer order(nu, nu, 2), bc(nu, 2, 2), e(nu, 2, 2)
      common /cstak/ ds
      double precision ds(500)
      integer ice, istkgt, max0, i, j, l
      integer ipps, inow, imaord, inoold, is(1000)
      real rs(1000)
      logical allero, ls(1000)
      complex cs(500)
      double precision ws(500)
      integer temp, temp1
      equivalence (ds(1), cs(1), ws(1), rs(1), is(1), ls(1))
c To determine which ODE should use which boundary condition.
c Mnemonic - Double precision Galerkin's method for Linear Systems,
c            Boundary condition Placement.
c Scratch Space Allocated -
c       S(DglsBP) <= Nu*(4*Nu+15)
c Integer words.
c Define Node = Is(inow) -> Nodei
c Check the input for errors.
      if (nu .lt. 1) call seterr(16hDglsBP - Nu.lt.1, 16, 1, 2)
      do  3 l = 1, 2
         do  2 i = 1, nu
c Is Order(i,.,l) = (-1,...,-1)?
            allero = .true.
            do  1 j = 1, nu
               allero = allero .and. order(i, j, l) .eq. (-1)
               if (order(i, j, l) .lt. (-1) .or. order(i, j, l) .gt. 2) 
     1            call seterr(
     2            41hDglsBP - Order(i,j,l) not one of -1,0,1,2, 41, 2, 2
     3            )
   1           continue
            if (bc(i, 1, l) .ne. (-2) .and. bc(i, 1, l) .ne. 0 .or. bc(i
     1         , 2, l) .ne. (-2) .and. bc(i, 2, l) .ne. 1) call seterr(
     2         36hDglsBP - BC(i,.,l) not one of -2,0,1, 36, 3, 2)
            if (allero) call seterr(
     1         33hDglsBP - Order(i,.,l)=(-1,...,-1), 33, 4, 2)
   2        continue
   3     continue
      call enter(1)
c Complement of E.
      ice = istkgt(nu, 2)
c Maxord(i,l) = Max over j=1,...,Nu Order(i,j,l).
      imaord = istkgt(2*nu, 2)
      call seti(2*nu, -1, is(imaord))
      do  6 l = 1, 2
         do  5 i = 1, nu
            do  4 j = 1, nu
               temp1 = imaord+i-1+(l-1)*nu
               temp = imaord+i-1+(l-1)*nu
               is(temp1) = max0(is(temp), order(i, j, l))
   4           continue
   5        continue
   6     continue
      i = 0
      ipps = 1
   7  if (i .ge. 4*nu .and. ipps .eq. 1) goto  15
         goto  12
c Make a node.
   8        inoold = inow
            i = i+1
            inow = istkgt(nu+3, 2)
            is(inow) = inoold
c       Get the candidates for E(i).
            call d6lsbp(i, nu, order, bc, e, is(imaord), is(ice), is(
     1         inow+3), is(inow+1))
            is(inow+2) = 0
            ipps = 0
            goto  13
            goto  13
c Searching a node.
   9        is(inow+2) = is(inow+2)+1
            if (is(inow+2) .le. is(inow+1)) goto 10
               ipps = -1
c Back-up.
               goto  7
  10        temp = inow+2+is(inow+2)-1
            e(i, 1, 1) = is(temp+1)
            ipps = 1
            goto  13
            goto  13
c Backing up a Node.
  11        inow = is(inow)
            call istkrl(1)
            i = i-1
            ipps = 0
            goto  13
            goto  13
  12        temp = ipps+2
            if (temp .gt. 0 .and. temp .le. 3) goto ( 11,  9,  8), temp
c End Switch.
  13     if (i .ne. 0) goto 14
            call seterr(37hDglsBP - Improper Boundary Conditions, 37, 5,
     1         1)
            goto  15
  14     continue
         goto  7
c End While.
  15  call leave
      return
      end
      subroutine d6lsbp(i, nu, order, bc, e, maxord, ce, r, n)
      integer nu
      integer i, order(nu, nu, 2), bc(1), e(1), maxord(nu, 2), ce(nu)
      integer r(nu), n
      integer mod, max0, j, l, nbcs, dm
      integer ii, lr
      integer temp
c BC(Nu,2,2),E(Nu,2,2),
c E(i-1),R(N).
      if (bc(i) .ge. 0) goto 1
         n = 1
         r(n) = 0
         return
c LR = 1 for left, LR = 2 for right.
   1  lr = (i-1)/(2*nu)+1
c DM = 1 for Dirichlet, DM = 2 for Mixed boundary conditions.
      dm = mod((i-1)/nu, 2)+1
      ii = mod(i, nu)
      if (ii .eq. 0) ii = nu
c B(i) = B(ii,DM,LR).
      n = 0
      do  2 j = 1, nu
         ce(j) = j
   2     continue
c CE = Complement of E.
      if (i .gt. 2*nu) goto 7
         j = 1
            goto  4
   3        j = j+1
   4        if (j .ge. i) goto  6
            if (bc(j) .lt. 0) goto 5
               temp = e(j)
               ce(temp) = 0
   5        continue
            goto  3
   6     continue
         goto  12
   7     j = 2*nu+1
            goto  9
   8        j = j+1
   9        if (j .ge. i) goto  11
            if (bc(j) .lt. 0) goto 10
               temp = e(j)
               ce(temp) = 0
  10        continue
            goto  8
  11     continue
  12  do  18 j = 1, nu
         if (ce(j) .eq. 0) goto  18
         nbcs = 0
         l = 1
            goto  14
  13        l = l+1
  14        if (l .ge. i) goto  15
            if (e(l) .eq. j .and. bc(l) .ge. 0) nbcs = nbcs+1
            goto  13
  15     if ((dm .ne. 1 .or. maxord(j, lr) .le. bc(i)) .and. (dm .ne. 2
     1       .or. order(j, ii, lr) .le. bc(i))) goto 17
            if (nbcs .ge. max0(maxord(j, 1), maxord(j, 2))) goto 16
               n = n+1
               r(n) = j
  16        continue
  17     continue
  18     continue
      return
      end