function [ adiag, aleft, arite, rhs ] = assemble_picard ( u, h, indx, n, nl, ...
  node, nquad, nu, problem, ul, ur, xn, xquad )

%% ASSEMBLE_PICARD assembles the Picard lienar system.
%
%  Discussion:
%
%    The equation we are trying to solve has the form:
%
%      -d/dx ( p(x) du/dx ) + q(x) * u + u * u' = f(x)
%
%    For the Picard iteration, we need to modify the nonlinear term u * u'
%    so that it is linear in the unknown u, and any other factors of u are
%    lagged.  One way to do this gives us the following equation:
%
%      -d/dx ( p(x) du/dx ) + q(x) * u + u * uold' = f(x)
%
%    where uold is the previous iterate.
%
%    Now we can formulate this system as a (linear) finite element problem
%
%      A * u = rhs
%
%    to be solved for the new approximate solution u.
%
%  Modified:
%
%    28 April 2007
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Parameters:
%
%    Input, real U(NU), the approximate solution from the previous 
%    iteration.
%
%    Input, real H(N), the length of the subintervals.  
%
%    Input, integer INDX(1:N+1).
%    For a node I, INDX(I) is the index of the unknown
%    associated with node I.
%    If INDX(I) is equal to -1, then no unknown is associated
%    with the node, because a boundary condition fixing the
%    value of U has been applied at the node instead.
%    Unknowns are numbered beginning with 1.
%    If IBC is 2 or 4, then there is an unknown value of U
%    at node 0, which will be unknown number 1.  Otherwise,
%    unknown number 1 will be associated with node 1.
%    If IBC is 1 or 4, then there is an unknown value of U
%    at node N, which will be unknown N or N+1,
%    depending on whether there was an unknown at node 0.
%
%    Input, integer N.
%    The number of subintervals into which the interval
%    [XL,XR] is broken.
%
%    Input, integer NL.
%    The number of basis functions used in a single
%    subinterval.  (NL-1) is the degree of the polynomials
%    used.  For this code, NL is fixed at 2, meaning that
%    piecewise linear functions are used as the basis.
%
%    Input, integer NODE(NL,N).
%    For each subinterval I:
%    NODE(1,I) is the number of the left node, and
%    NODE(2,I) is the number of the right node.
%
%    Input, integer NQUAD.
%    The number of quadrature points used in a subinterval.
%    This code uses NQUAD = 1.
%
%    Input, integer NU.
%    NU is the number of unknowns in the linear system.
%    Depending on the value of IBC, there will be N-1,
%    N, or N+1 unknown values, which are the coefficients
%    of basis functions.
%
%    Input, integer PROBLEM, indicates which problem to be solved.
%    * 1, u = x, p=1, q=0, f=x, u(0)=0, u'(1)=1.
%    * 2, u = 2*(1-cos(0.5*pi*x))/pi, p=1, q=0,
%    f = -0.5*pi*cos(0.5*pi*x) + 2*sin(0.5*pi*x)*(1-cos(0.5*pi*x)/pi
%    u(0) = 0, u'(1)=1.
%
%    Input, real UL.
%    If IBC is 1 or 3, UL is the value that U is required
%    to have at X = XL.
%    If IBC is 2 or 4, UL is the value that U' is required
%    to have at X = XL.
%
%    Input, real UR.
%    If IBC is 2 or 3, UR is the value that U is required
%    to have at X = XR.
%    If IBC is 1 or 4, UR is the value that U' is required
%    to have at X = XR.
%
%    Input, real XN(1:N+1).
%    XN(I) is the location of the I-th node.  XN(1) is XL,
%    and XN(N+1) is XR.
%
%    Input, real XQUAD(N)
%    XQUAD(I) is the location of the single quadrature point
%    in interval I.
%
%    Output, real ADIAG(NU).
%    ADIAG(I) is the "diagonal" coefficient of the I-th
%    equation in the linear system.  That is, ADIAG(I) is
%    the coefficient of the I-th unknown in the I-th equation.
%
%    Output, real ALEFT(NU).
%    ALEFT(I) is the "left hand" coefficient of the I-th
%    equation in the linear system.  That is, ALEFT(I) is the
%    coefficient of the (I-1)-th unknown in the I-th equation.
%    There is no value in ALEFT(1), since the first equation
%    does not refer to a "0-th" unknown.
%
%    Output, real ARITE(NU).
%    ARITE(I) is the "right hand" coefficient of the I-th
%    equation in the linear system.  ARITE(I) is the coefficient
%    of the (I+1)-th unknown in the I-th equation.  There is
%    no value in ARITE(NU) because the NU-th equation does not
%    refer to an "NU+1"-th unknown.
%
%    Output, real RHS(NU), the right hand side of the linear
%    equations.
%
  rhs(1:nu) = 0.0;
  adiag(1:nu) = 0.0;
  aleft(1:nu) = 0.0;
  arite(1:nu) = 0.0;
%
%  For element IE...
%
  for ie = 1 : n

    he = h(ie);
    xleft = xn(node(1,ie)+1);
    xrite = xn(node(2,ie)+1);
%
%  For quadrature point IQ...
%
    for iq = 1 : nquad

      xqe = xquad(ie);
%
%  Compute value of U for previous solution.
%
      total = 0.0;

      for il = 1 : nl

        ig = node(il,ie);
        iu = indx(ig+1);

        if ( iu <= 0 )

          if ( il == 1 )
            total = total + ul;
          else
            total = total + ur;
          end

        else
          total = total + u(iu);
        end

      end

      uold = total / nl;
%
%  Compute value of U' for previous solution.
%
      jl = node(1,ie);
      jr = node(2,ie);
      iul = indx(jl+1);
      iur = indx(jr+1);

      if ( iul <= 0 )
        uoldx = ( u(iur) - ul ) / he;
      elseif ( iur <= 0 )
        uoldx = ( ur - u(iul) ) / he;
      else
        uoldx = ( u(iur) - u(iul) ) / he;
      end
%
%  For basis function IL...
%
      for il = 1 : nl

        ig = node(il,ie);
        iu = indx(ig+1);

        if ( 0 < iu )

          [ phii, phiix ] = phi ( il, xqe, xleft, xrite );

          rhs(iu) = rhs(iu) + he * phii * ff ( xqe, problem );
%
%  Handle boundary conditions that prescribe the value of U'.
%
          if ( ig == 0 )

            x = 0.0;
            rhs(iu) = rhs(iu) - pp ( x, problem ) * ul;

          elseif ( ig == n )

            x = 1.0;
            rhs(iu) = rhs(iu) + pp ( x, problem ) * ur;

          end
%
%  For basis function JL...
%
          for jl = 1 : nl

            jg = node(jl,ie);
            ju = indx(jg+1);

            [ phij, phijx ] = phi ( jl, xqe, xleft, xrite );

            aij = he * ( pp ( xqe, problem ) * phiix * phijx ...
                  + qq ( xqe, problem ) * phii * phij ...
                  + uold * phii * phijx );

            if ( ju <= 0 )

              if ( jg == 0 )
                rhs(iu) = rhs(iu) - aij * ul;
              elseif ( jg == n )
                rhs(iu) = rhs(iu) - aij * ur;
              end

            elseif ( iu == ju )
              adiag(iu) = adiag(iu) + aij;
            elseif ( ju < iu )
              aleft(iu) = aleft(iu) + aij;
            else
              arite(iu) = arite(iu) + aij;
            end

          end

        end

      end

    end

  end