function ortho ( a, alpha, beta, np, problem, quad_num, quad_w, quad_x )

%% ORTHO tests the basis functions for orthogonality.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 November 2006
%
%  Author:
%
%    Original FORTRAN77 version by Max Gunzburger, Teresa Hodge
%    MATLAB version by John Burkardt
%
%  Parameters:
%
%    Input, real A(1:NP+1), the squares of the norms of the 
%    basis functions.
%
%    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, integer NP.
%    The highest degree polynomial to use.
%
%    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.
%

%
%  Zero out the B array, so we can start summing up the dot products.
%
  b(1:np+1,1:np+1) = 0.0;
%
%  Approximate the integral of the product of basis function
%  I and basis function J over the interval [-1,1].
%
%  We expect to get zero, except when I and J are equal,
%  when we should get A(I).
%
  for iq = 1 : quad_num
    x = quad_x(iq);
    for i = 0 : np
      [ phii, phiix ] = phi ( alpha, beta, i, np, x );
      for j = 0 : np

        [ phij, phijx ] = phi ( alpha, beta, j, np, x );

        bij = pp ( x, problem ) * phiix * phijx ...
            + qq ( x, problem ) * phii * phij;

        b(i+1,j+1) = b(i+1,j+1) + bij * quad_w(iq);

      end
    end
  end
%
%  Print out the results of the test.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'Basis function orthogonality test:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '   i   j     b(i,j)/a(i)\n' );
  fprintf ( 1, '\n' );
  for i = 0 : np
    fprintf ( 1, '\n' );
    for j = 0 : np
      fprintf ( 1, '  %4d  %4d  %12f\n', i, j, b(i+1,j+1) / a(i+1) );
    end
  end

  return
end