function [ a, alpha, beta ] = alpbet ( np, problem, quad_num, quad_w, quad_x )

%% ALPBET calculates the coefficients in the recurrence relationship.
%
%  Discussion:
%
%    ALPHA and BETA are the coefficients in the three
%    term recurrence relation for the orthogonal basis functions
%    on [-1,1].
%
%    The routine also calculates A, the square of the norm of each basis
%    function.
%
%  Modified:
%
%    03 November 2006
%
%  Author:
%
%    Max Gunzburger
%    Teresa Hodge
%
%    MATLAB translation by John Burkardt
%
%  Parameters:
%
%    Output, real A(1:NP+1), the squares of the norms of the 
%    basis functions.
%
%    Output, real ALPHA(NP).
%    ALPHA(I) contains one of the coefficients of a recurrence
%    relationship that defines the basis functions.
%
%    Output, 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.
%
  ss = 0.0;
  su = 0.0;

  for iq = 1 : quad_num

    x = quad_x(iq);

    s = 4.0 * pp ( x, problem ) * x * x ...
      + qq ( x, problem ) * ( 1.0 - x * x )^2;

    u = 2.0 * pp ( x, problem ) * x * ( 3.0 * x * x - 1.0 ) ...
      + x * qq ( x, problem ) * ( 1.0 - x * x )^2;

    ss = ss + s * quad_w(iq);
    su = su + u * quad_w(iq);

  end

  beta(1) = 0.0;
  a(1) = ss;
  alpha(1) = su / ss;

  for i = 2 : np+1

    ss = 0.0;
    su = 0.0;
    sv = 0.0;

    for iq = 1 : quad_num

      x = quad_x(iq);
      q = 1.0;
      qm1 = 0.0;
      qx = 0.0;
      qm1x = 0.0;

      for k = 1 : i - 1
        qm2 = qm1;
        qm1 = q;
        qm2x = qm1x;
        qm1x = qx;
        q = ( x - alpha(k) ) * qm1 - beta(k) * qm2;
        qx = qm1 + ( x - alpha(k) ) * qm1x - beta(k) * qm2x;
      end

      t = 1.0 - x * x;

      s = pp ( x, problem ) * ( t * qx - 2.0 * x * q )^2 ...
        + qq ( x, problem ) * ( t * q )^2;

      u = pp ( x, problem ) ...
        * ( x * t * qx + ( 1.0 - 3.0 * x * x ) * q ) ...
        * ( t * qx - 2.0 * x * q ) + x * qq ( x, problem ) * ( t * q )^2;

      v = pp ( x, problem ) ...
        * ( x * t * qx + ( 1.0 - 3.0 * x * x ) * q ) ...
        * ( t * qm1x - 2.0 * x * qm1 ) ...
        + x * qq ( x, problem ) * t * t * q * qm1;

      ss = ss + s * quad_w(iq);
      su = su + u * quad_w(iq);
      sv = sv + v * quad_w(iq);

    end

    a(i) = ss;

    if ( i <= np )
      alpha(i) = su / ss;
      beta(i) = sv / a(i-1);
    end

  end