function f = sol ( a, alpha, beta, np, problem, quad_num, quad_w, quad_x )

%% SOL solves a linear system for the finite element coefficients.
%
%  Modified:
%
%    03 November 2006
%
%  Author:
%
%    Max Gunzburger
%    Teresa Hodge
%
%    MATLAB translation 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.
%
%    Output, 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.
%
  f(1:np+1) = 0.0;

  for iq = 1 : quad_num
    x = quad_x(iq);
    t = ff ( x, problem ) * quad_w(iq);
    for i = 0 : np
      [ phii, phiix ] = phi ( alpha, beta, i, np, x );
      f(i+1) = f(i+1) + phii * t;
    end
  end

  f(1:np+1) = f(1:np+1) ./ a(1:np+1);