% Gauss Quadrature for Example 11.9
clearvars
clc

g = @(z) z.*exp(z);
a = 0;
b = 3;

Iexact = integral(g,a,b);

% Change of variables
f = @(x) g(((b-a).*x + b + a)/2).*(b-a)/2;

% Gauss Quadrature with transformed function
I1 = gauss_quad(f,1);
I2 = gauss_quad(f,2);
I3 = gauss_quad(f,3);
I4 = gauss_quad(f,4);

Er1 = (Iexact - I1)/Iexact;
Er2 = (Iexact - I2)/Iexact;
Er3 = (Iexact - I3)/Iexact;
Er4 = (Iexact - I4)/Iexact;

disp('GAUSS QUADRATURE')
disp('------------------------------------------------------------')
disp(['With 1 integration point,  I1 = ',num2str(I1),' and Er1 = ',num2str(Er1)])
disp(['With 2 integration points, I2 = ',num2str(I2),' and Er2 = ',num2str(Er2)])
disp(['With 3 integration points, I3 = ',num2str(I3),' and Er3 = ',num2str(Er3)])
disp(['With 4 integration points, I4 = ',num2str(I4),' and Er4 = ',num2str(Er4)])

xplot = linspace(a,b);
yplot = g(xplot);

figure
hold on
fplot(g,[a,b],'LineWidth',2)
area(xplot,yplot,'EdgeColor',[0.8500, 0.3250, 0.0980],'FaceAlpha',0.2)
xlim([a,b])
plot([0,0],ylim,'-k')
plot(xlim,[0,0],'-k')
xlabel('x')
ylabel('f(x)')
title('f(x) = xe^x')