Implement GaussLegendre Integration in MATLAB Create a funti

Implement Gauss-Legendre Integration in MATLAB. Create a funtion called \'integrateGL\'. The function should have four input arguments: a function, lower bound of integration interval, upper bound of intergration interval, and mesh points N. The mesh points are of the set [2, 4, 6, 8, 10, 12, 16, 24].
.dat files will be imported into MATLAB and used as the associated weights for each number of mesh points. (The contents of the .dat files, for example, are as follows)
GAUSS - 2
0.577350269189626e0 1e0
-0.577350269189626e0 1e0

GAUSS - 4
0.861136311594053e0 0.347854845137454e0
0.339981043584856e0 0.652145154862546e0
-0.339981043584856e0 0.652145154862546e0
-0.861136311594053e0 0.347854845137454e0

etcetera

Solution

CODE:

function [x,w]=integrateGL(N,a,b)

% Computes the Legendre-Gauss nodes and weights on an interval
% [a,b] with truncation order N
%
% For a continuous function f(x) which is defined on [a,b]
% if we can evaluate at any x in [a,b]. Simply evaluate it at all of
% the values contained in the x vector to obtain a vector f.

% Then compute the definite integral using sum(f.*w);

clc;
N=N-1;
N1=N+1; N2=N+2;

xu=linspace(-1,1,N1)\';

% Initial guess
y=cos((2*(0:N)\'+1)*pi/(2*N+2))+(0.27/N1)*sin(pi*xu*N/N2);

% Create an Legendre-Gauss Vandermonde Matrix
L=zeros(N1,N2);

% Derivative of LGVM
Lp=zeros(N1,N2);

% Compute the zeros of the N+1 Legendre Polynomial
% using the Newton-Raphson method and the recursion relation.
y0=2;

% Iterate till the new points are uniformly within epsilon of old points
while max(abs(y-y0))>eps
Lp(:,1)=0;
L(:,1)=1;
Lp(:,2)=1;
L(:,2)=y;
  
for i=2:N1
L(:,i+1)=( (2*i-1)*y.*L(:,i)-(i-1)*L(:,i-1) )/i;
end

Lp=(N2)*( L(:,N1)-y.*L(:,N2) )./(1-y.^2);   
y0=y;
y=y0-L(:,N2)./Lp;
end

% Linear map from [-1,1] to [a,b]
x=(a*(1-y)+b*(1+y))/2;

% Compute the weights
w=(b-a)./((1-y.^2).*Lp.^2)*(N2/N1)^2;

end