function [ a, f ] = assemble_backward_euler ( node_num, node_xy, ...
  element_order, element_num, element_node, quad_num, ib, time, ...
  time_step_size, u_old, a, f )

%% ASSEMBLE_BACKWARD_EULER adjusts the system for the backward Euler term.
%
%  Discussion:
%
%    The input linear system
%
%      A * U = F
%
%    is appropriate for the equation
%
%      -Uxx - Uyy - K * U = RHS
%
%    We need to modify the matrix A and the right hand side F to
%    account for the approximation of the time derivative in
%
%      Ut - Uxx - Uyy - K * U = RHS
%
%    by the backward Euler approximation:
%
%      Ut approximately equal to ( U - Uold ) / dT
%
%  Modified:
%
%    22 July 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NODE_NUM, the number of nodes.
%
%    Input, real NODE_XY(2,NODE_NUM), the coordinates of nodes.
%
%    Input, integer ELEMENT_ORDER, the number of nodes used to form one element.
%
%    Input, integer ELEMENT_NUM, the number of elements.
%
%    Input, integer ELEMENT_NODE(ELEMENT_ORDER,ELEMENT_NUM);
%    ELEMENT_NODE(I,J) is the global index of local node I in element J.
%
%    Input, integer QUAD_NUM, the number of quadrature points used in assembly.
%
%    Input, integer IB, the half-bandwidth of the matrix.
%
%    Input, real TIME, the current time.
%
%    Input, real TIME_STEP_SIZE, the size of the time step.
%
%    Input, real U_OLD(NODE_NUM), the finite element
%    coefficients for the solution at the previous time.
%
%    Input, real A(3*IB+1,NODE_NUM), the NODE_NUM
%    by NODE_NUM coefficient matrix, stored in a compressed format.
%
%    Input, real F(NODE_NUM), the right hand side.
%
%    Output, real A(3*IB+1,NODE_NUM), the updated matrix.
%
%    Output, real F(NODE_NUM), the updated right hand side.
%

%
%  Get the quadrature rule weights and nodes.
%
  [ quad_w, quad_xy ] = quad_rule ( quad_num );

  for element = 1 : element_num
%
%  Make two copies of the triangle.
%
    t3(1:2,1:3) = node_xy(1:2,element_node(1:3,element));
    t6(1:2,1:6) = node_xy(1:2,element_node(1:6,element));
%
%  Map the quadrature points QUAD_XY to points PHYS_XY in the physical triangle.
%
    phys_xy = reference_to_physical_t3 ( t3, quad_num, quad_xy );

    area = abs ( triangle_area_2d ( t3 ) );

    w(1:quad_num) = area * quad_w(1:quad_num);

    for quad = 1 : quad_num

      for test = 1 : element_order

        node = element_node(test,element);

        [ bi, dbidx, dbidy ] = basis_11_t6 ( t6, test, phys_xy(1:2,quad) );
%
%  Carry the U_OLD term to the right hand side.
%
        f(node) = f(node) + w(quad) * bi * u_old(node) / time_step_size;
%
%  Modify the diagonal entries of A.
%
        for basis = 1 : element_order

          j = element_node(basis,element);

          [ bj, dbjdx, dbjdy ] = basis_11_t6 ( t6, basis, phys_xy(1:2,quad) );

          a(node-j+2*ib+1,j) = a(node-j+2*ib+1,j) ...
            + w(quad) * bi * bj / time_step_size;

        end
      end

    end

  end