function l2_error = compute_frac_l2_error ( beta, coeff_end, n, x, w1_h , indicator)     

%*****************************************************************************80
%
%% COMPUTE_FRAC_L2_ERROR estimates the L2 error norm 
%
%
%  Parameters:
%
%    Input, integer N, the number of nodes.
%
%    Input, real X(N), the mesh points.
%
%    Input, real U(N), the finite element coefficients.
%
%    Output, real L2_ERROR, the estimated L2 norm of the error.
%
  l2_error = 0.0;
% Mesh size:
  h = 1/(n-1);
  
% True Solution of u:
  true_u = @(x) x.^2.*(1-x).^2; 
  
  if indicator == 1       % True Solution of DerBeta_w2: (w2 = x) 
      true_DerBeta_w2 = @(x) 1/gamma(2-beta)*x.^(1-beta);    
  else                    % True Solution of DerBeta_w2: (w2 = (x^2+2*x)/3))
      true_DerBeta_w2 = @(x) 2/(3*gamma(3-beta))*(x.^(2-beta) + (2-beta)*x.^(1-beta));    
  end
%
%  Quadrature definitions.
%
  quad_num = 2;
  abscissa(1) = -0.577350269189625764509148780502;
  abscissa(2) = +0.577350269189625764509148780502;
  weight(1) = 1.0;
  weight(2) = 1.0;
%
%  Integrate over each interval.
%
  for i = 1 : n - 1

    xl = x(i);
    xr = x(i+1);
    w1l = w1_h(i);
    w1r = w1_h(i+1);

    for q = 1 : quad_num

      xq = ( ( 1.0 - abscissa(q) ) * xl   ...
           + ( 1.0 + abscissa(q) ) * xr ) ...
           /   2.0;

      wq = weight(q) * ( xr - xl ) / 2.0;
%
%  Compute DerBeta_w1 at xq.
%
     if i == 1
         DerBeta_w1q = w1r/(gamma(2-beta)*h)*xq^(1-beta);
     else
        temp0 = (xq-x(2:i)+h).^(1-beta) - (xq-x(2:i)).^(1-beta);
        temp1 = temp0.*w1_h(2:i);
        temp2 = temp0.*w1_h(1:i-1);
        temp = sum(temp1) - sum(temp2);
        DerBeta_w1q = 1/(gamma(2-beta)*h)*(temp + (w1r-w1l)*(xq-xl)^(1-beta));          
     end
     
      DerBeta_w2q = true_DerBeta_w2 ( xq );    
      uq = DerBeta_w1q + coeff_end*DerBeta_w2q ;
      eq = true_u(xq);

      l2_error = l2_error + wq * ( uq - eq )^2;
     end
  end

  l2_error = sqrt ( l2_error );

  return
end
