echo off;
% NUMERICAL METHODS: MATLAB Programs, (c) John H. Mathews 1994
% To accompany the text:
% NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992
% Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A.
% This free software is complements of the author.

% Algorithm 7.6 (Gauss-Legendre Quadrature).
% Section 7.5, Gauss-Legendre Integration, Page 397
echo on; clc; format long; hold off; clear

% This program implements Gaussian quadrature.

% The function f(x) is integrated over the interval [a,b].

% Define and store the function f(x) in the M-file f.m
% function y = f(x)
% y = 13.*(x - x.^2).*exp(-3.*x./2);

delete f.m
diary f.m; disp('function y = f(x)');...
disp('y = 13.*(x - x.^2).*exp(-3.*x./2);');...
diary off;

f(0); % Remark. f.m gauss.m gauraw.m gauaw.mat are used for Algorithm 7.6

pause % Press any key to see the graph y = f(x).

% The abscissas and weights are loaded from the file gauaw.mat

% gauraw.m contains the instructions to read this file.

clc; clg;
a = 0;
b = 4;
h = (b-a)/200;
X1 = a:h:b;
Y1 = f(X1);
c = -1.5;
d = 2;
axis([a b c d]);...
plot(X1,Y1,'g');...
hold on;...
plot([a b],[0 0],'b',[0 0],[c d],'b');...
xlabel('x');...
ylabel('y');...
title('The curve y = f(x) over [a,b].');...
grid;...
axis;...
hold off;...
shg; pause % Press any key to continue.

clc;

% Be patient for a moment.

[A,W] = gauraw; % Loading the abscissas and weights.

pause % Press any key to continue.

clc;

% Gaussian quadrature is applied over the interval [a,b].

% Place the endpoints of [a,b] in a and b.

% Place the tolerance in toler.

a = 0;
b = 4;
toler = 0.00001;

[Q,N,X,quad,err] = gauss('f',a,b,A,W,toler);

pause % Press any key to continue.

clg;
Y = f(X);
Z = zeros(1,length(X));
c = -1.5;
d = 2;
axis([a b c d]);...
plot(X1,Y1,'g',X,Y,'or',X,Z,'+r');...
hold on;...
plot([a b],[0 0],'b',[0 0],[c d],'b');...
xlabel('x');...
ylabel('y');...
title('Gaussian quadrature.');...
grid;...
axis;...
hold off;...
shg; pause % Press any key to continue.

values = [N;Q];
Mx1 = 'The Gaussian quadrature approximations: ';
Mx2 = ' No. of nodes quadrature value';
Mx3 = 'The approximate value of the integral is:';
Mx4 = ' quadrature value +- error bound';
clc,echo off,diary output,...
disp(Mx1),disp(''),disp(Mx2),disp(values'),...
disp(Mx3),disp(Mx4),disp([quad err]),,diary off,echo on