% EXAMPLE2_33 Computation for Example 2.33.
%
%    At the bottom of the routine quadgp.m, the weight
%    function wf has to be specified as follows:
%           y=x.^x.*(1-x).^(1-x);
%
clear mc uv AB
global mc uv AB
eps0=.5e-14;
mc=3; exact=.617826964729011473;
AB=[[0 .00115];[.00115 .0844];[.0844 .5]];
n=10; int=0;
while(abs(int-exact)>eps0)
  sm=0;
  uv=fejer(n);
  for i=1:3
    xw=quadgp(n,i);
    sm=sm+sum(xw(:,2));
  end
  int=2*sm;
  n=n+1;
end 
n
int
err=abs(int-exact)
mc=1; AB=[0 .5]; i=1;
uv=fejer(108);
xw=quadgp(108,i);
sm=sum(xw(:,2));
int=2*sm
err=abs(int-exact)
uv=fejer(324);
xw=quadgp(324,i);
sm=sum(xw(:,2));
int=2*sm
err=abs(int-exact)