function [ n, xn, status ] = subdiv ( eta, kount, n, nmax, tol, xn )

%% SUBDIV decides which intervals should be subdivided.
%
%  Modified:
%
%    06 November 2006
%
%  Parameters:
%
%    Input, real ETA(N).
%    ETA(I) is the error estimate for interval I.  It is computed
%    as the sum of two quantities, one associated with the left
%    and one with the right node of the interval.
%
%    Input, integer KOUNT, the number of adaptive steps that have been taken.
%
%    Input, integer N
%    The number of subintervals into which the interval
%    [XL,XR] is broken.  
%
%    Input, integer NMAX, the maximum number of unknowns that can be handled.
%
%    Input, real TOL.
%    A tolerance that is used to determine whether the estimated
%    error in an interval is so large that it should be subdivided
%    and the problem solved again.
%
%    Input, real XN(0:N).
%    XN(I) is the location of the I-th node.  XN(0) is XL,
%    and XN(N) is XR.
%
%    Output, integer N, the updated number of nodes.
%
%    Output, real XN(0:N), the updated nodes.
%
%    Output, integer STATUS, reports status of subdivision.
%    0, a new subdivision was carried out.
%    1, no more subdivisions are needed.
%    -1, no more subdivisions can be carried out.
%
%  Local Parameters:
%
%    Local, integer JADD(N).
%    JADD(I) is 1 if the error estimates show that interval I
%    should be subdivided.
%
  status = 0;
%
%  Add up the ETA's, and get their average.
%
  ave = sum ( eta(1:n) ) / n;
%
%  Look for intervals whose ETA value is relatively large,
%  and note in JADD that these intervals should be subdivided.
%
  k = 0;
  temp = max ( 1.2 * ave + 0.00001, tol^2 / n );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'Tolerance = %f\n', temp );
  fprintf ( 1, '\n' );

  for j = 1 : n

    if ( temp < eta(j) )
      k = k + 1;
      jadd(j) = 1;
      fprintf ( 1, 'Subdivide interval %d\n', j );
    else
      jadd(j) = 0;
    end

  end
%
%  If no subdivisions needed, we're done.
%
  if ( k <= 0 )
    fprintf ( 1, 'Success on step %d\n', kount );
    status = 1;
    return
  end
%
%  See if we're about to go over our limit.
%
  if ( nmax < n + k )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'The iterations did not reach their goal.\n' );
    fprintf ( 1, 'The next value of N is %d\n', n + k );
    fprintf ( 1, 'which exceeds NMAX = %d\n', nmax );
    status = -1;
    return
  end
%
%  Insert new nodes where needed.
%
  k = 0;
  xt(1) = xn(1);
  for j = 1 : n

    if ( 0 < jadd(j) )
      xt(j+k+1) = 0.5 * ( xn(j+1) + xn(j) );
      k = k + 1;
    end

    xt(j+k+1) = xn(j+1);

  end
%
%  Update the value of N, and copy the new nodes into XN.
%
  n = n + k;

  xn(1:n+1) = xt(1:n+1);