parameter (nx=7, nq=5)
      implicit double precision (a-g,o-z)
      double precision f(nx), df(nx,nx), work(nx), 
     c        xa(nx), ax(nq,(nx+1)/2), wx(nq,(nx+1)/2),c0(nx), c1(nx)
      integer iwork(nx), nxp((nx+1)/2), ibc(nx)
c
      na   = nx
      nxm1 = nx-1
      ne   = (nx+1)/2-1
      nbc  = 1
c
c     Open input and output files.
c
      open(unit=8, file='junk.dat')
c
c     Read in parameters, grid, boundary conditions, and if specified, 
c     the initial guess for consumption.
c
      xa(1) = 0.0d0
      xa(2) = 0.5d0
      xa(3) = 1.0d0
      xa(4) = 1.5d0
      xa(5) = 3.0d0
      xa(6) = 4.5d0
      xa(7) = 6.0d0

      c0(1) = 1.0
      c0(2) = 0.6065
      c0(3) = 0.3679
      c0(4) = 0.2231
      c0(5) = 0.0498
      c0(6) = 0.0111
      c0(7) = 0.0025
      
      ibc(1) = 1
      ibc(2) = 0
      ibc(3) = 0
      ibc(4) = 0
      ibc(5) = 0
      ibc(6) = 0
      ibc(7) = 0

      do 80 i=1,ne
        nxp(i) = 3
 80   continue
      crit  = 0.0000001d0
      maxit = 20
c
c     Compute quadrature abscissas and weights.
c
      do 110 i=1,ne
        call qgausl (nxp(i), xa(2*i-1), xa(2*i+1), ax(1,i), wx(1,i))
 110  continue
c
c     The Newton-Raphson iteration, i.e., find x such that f(x)=0 where 
c     x are coefficients of the finite element approximation.
c
      do 130 it=1,maxit
c
c       Initialize f and its derivative df.
c
        do 140 i=1,na
          f(i) = 0.d0
          do 140 j=1,na
            df(i,j) = 0.d0
 140    continue
c
c       Construct f(c0) = int R(x,j;c0) N_n(x) dx and its derivative on
c       each element n.
c
          do 160 n=1,ne
c
c           Element n is the interval [xa(n),xa(n+1)] 
c
            x1   = xa(2*n-1) 
            dx   = xa(2*n+1)-x1 
            a    = (xa(2*n)-x1)/dx 
            a1   = 1.d0-a
c
c           Compute f at all quadrature points on the element.
c
            do 170 i=1,nxp(n)
              x     = ax(i,n)
              u     = wx(i,n)
              xi    = (x-x1)/dx
              bs1   = (1.d0-xi)*(1.d0-xi/a)
              bs2   = xi/(a*a1)*(1.d0-xi)
              bs3   = -a*xi/a1*(1.d0-xi/a)
              dxi   = 1.d0/dx
              dbs1  = (-a-1.d0+2.d0*xi)*dxi/a
              dbs2  = (1.d0-2.d0*xi)*dxi/a/a1
              dbs3  = (-a+2.d0*xi)*dxi/a1
c
c             Compute function and residual using finite element approximation.
c
              c     = c0(2*n-1)*bs1 +c0(2*n)*bs2 +c0(2*n+1)*bs3
              dc    = c0(2*n-1)*dbs1+c0(2*n)*dbs2+c0(2*n+1)*dbs3
              res   = c+dc
c
c             Add the residual x quadrature weights x shape functions to f(.). 
c
              f(2*n-1)   = f(2*n-1)  + res*bs1*u
              f(2*n)     = f(2*n)    + res*bs2*u
              f(2*n+1)   = f(2*n+1)  + res*bs3*u
c
c             Add derivative of the residual x quadrature weights x shape
c             functions to df(.). 
c
              df(2*n-1,2*n-1) = df(2*n-1,2*n-1) + bs1*(bs1+dbs1)*u
              df(2*n-1,2*n)   = df(2*n-1,2*n)   + bs1*(bs2+dbs2)*u
              df(2*n-1,2*n+1) = df(2*n-1,2*n+1) + bs1*(bs3+dbs3)*u
              df(2*n,2*n-1)   = df(2*n,2*n-1)   + bs2*(bs1+dbs1)*u
              df(2*n,2*n)     = df(2*n,2*n)     + bs2*(bs2+dbs2)*u
              df(2*n,2*n+1)   = df(2*n,2*n+1)   + bs2*(bs3+dbs3)*u
              df(2*n+1,2*n-1) = df(2*n+1,2*n-1) + bs3*(bs1+dbs1)*u
              df(2*n+1,2*n)   = df(2*n+1,2*n)   + bs3*(bs2+dbs2)*u
              df(2*n+1,2*n+1) = df(2*n+1,2*n+1) + bs3*(bs3+dbs3)*u
 170        continue
 160      continue
c
c       Invert the derivative of f (i.e., df).
c
        ii=0
        do 230 i=1,na
          if (ibc(i).eq.0) then
            ii = ii+1
            f(ii) = f(i)
            do 240 j=1,na
              df(ii,j) = df(i,j)
 240        continue
            do 250 j=1,na
              df(j,ii) = df(j,i)
 250        continue
          endif
 230    continue
        call dgesv(na-nbc,1,df,na,iwork,f,na,info)
c
c       Issue warnings if necessary.
c
        if (info.lt.0) write(*,*) 'WARNING: illegal value for', 
     c    ' argument ',-info,' in DGESV'
        if (info.gt.0) write(*,*) 'WARNING: singular matrix',
     c    ' when solving df*x=f in SGM1DL'
c
c       Do Newton update.
c
        sum1 = 0.d0
        ii   = 0
        do 260 i=1,nx
          if (ibc(i).eq.0) then
            ii    = ii+1
            c1(i) = c0(i) - f(ii)
            sum1  = sum1 + f(ii)*f(ii)
          else
            c1(i) = c0(i)
          endif
 260    continue
        sum1  = dsqrt(sum1)/float(na-nbc)
c
c       Check to see if the solution is converged.   If so, stop.
c
        if (sum1.lt.crit) go to 999
        do 270 i=1,nx
            c0(i) =  c1(i) 
 270    continue
 130  continue
 999  continue
      do 320 i=1,nx
        write(8,1) xa(i), c1(i)
 320  continue
c
c     Format statements.
c
 1    format (9X,F8.3,1X,F8.6)
      stop
      end