function exact ( alpha, beta, f, np, nprint, problem, quad_num, ...
  quad_w, quad_x )

%% EXACT compares the computed and exact solutions.
%
%  Modified:
%
%    03 November 2006
%
%  Author:
%
%    Max Gunzburger
%    Teresa Hodge
%
%    MATLAB translation by John Burkardt
%
%  Parameters:
%
%    Input, real ALPHA(NP).
%    ALPHA(I) contains one of the coefficients of a recurrence
%    relationship that defines the basis functions.
%
%    Input, real BETA(NP).
%    BETA(I) contains one of the coefficients of a recurrence
%    relationship that defines the basis functions.
%
%    Input, real F(1:NP+1).
%    F contains the basis function coefficients that form the
%    representation of the solution U.  That is,
%      U(X)  =  SUM (I=0 to NP) F(I+1) * BASIS(I)(X)
%    where "BASIS(I)(X)" means the I-th basis function
%    evaluated at the point X.
%
%    Input, integer NP.
%    The highest degree polynomial to use.
%
%    Input, integer NPRINT.
%    The number of points at which the computed solution
%    should be printed out at the end of the computation.
%
%    Input, integer PROBLEM, indicates the problem being solved.
%    1, U=1-x**4, P=1, Q=1, F=1.0+12.0*x**2-x**4.
%    2, U=cos(0.5*pi*x), P=1, Q=0, F=0.25*pi*pi*cos(0.5*pi*x).
%
%    Input, integer QUAD_NUM, the order of the quadrature rule.
%
%    Input, real QUAD_W(QUAD_NUM), the quadrature weights.
%
%    Input, real QUAD_X(QUAD_NUM), the quadrature abscissas.
%
  nsub = 10;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'Comparison of computed and exact solutions:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    X        U computed    U exact     Difference\n' );
  fprintf ( 1, '\n' );

  for i = 0 : nprint
    x = ( 2 * i - nprint ) / nprint;
    ue = uex ( x, problem );
    up = 0.0;
    for j = 0 : np
      [ phii, phiix ] = phi ( alpha, beta, j, np, x );
      up = up + phii * f(j+1);
    end
    fprintf ( 1, '  %8f  %12f  %12f  %12f\n', x, up, ue, ue - up );
  end
%
%  Compute the big L2 error.
%
  big_l2 = 0.0;

  for i = 1 : nsub

    xl = ( 2 * i - nsub - 1 ) / nsub;
    xr = ( 2 * i - nsub     ) / nsub;

    for j = 1 : quad_num

      x = ( xl * ( 1.0 - quad_x(j) ) ...
          + xr * ( 1.0 + quad_x(j) ) ) / 2.0;

      up = 0.0;
      for k = 0 : np
        [ phii, phiix ] = phi ( alpha, beta, k, np, x );
        up = up + phii * f(k+1);
      end

      big_l2 = big_l2 + ( up - uex ( x, problem ) )^2 * quad_w(j) ...
        * ( xr - xl ) / 2.0;

    end

  end

  big_l2 = sqrt ( big_l2 );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'Big L2 error = %f\n', big_l2 );