M = 10;
xi = double(GLL_nodes(M));
a = 0;
b = 2*pi;
x = 0.5*(1-xi)*a+0.5*(1+xi)*b;
ya = 1;
yb = 1;
J = (b-a)/2;
D = DerivMatrix(xi,M)/J;
D2 = D*D;
L = D2+eye(M+1);
L_hat = L(2:end-1,2:end-1);
b_hat = zeros(M-1,1)-ya*L(2:end-1,1)-yb*L(2:end-1,end);
sol = L_hat\b_hat;
plot(x,[ya;sol;yb])
ylabel(‘y’)
xlabel(‘x’)
function D=DerivMatrix(x,n)
D=zeros(n+1,n+1);
for i=1:n+1
num(i)=1.0;
for k=1:n+1
if k ~= i
num(i)=num(i)*(x(i)-x(k));
end
end
end
for j=1:n+1
den(j)=1.0;
for k=1:n+1
if k ~= j
den(j)=den(j)*(x(j)-x(k));
end
end
end
for i=1:n+1
for j=1:n+1
if i ~= j
D(i,j)=num(i)/(den(j)*(x(i)-x(j)));
end
end
end
for i=1:n+1
for k=1:n+1
if i~=k
D(i,i)=D(i,i)+1./(x(i)-x(k));
end
end
end
end
function yint=LagrangePolynomial(x,y,n,xint)
yint=zeros(size(xint));
for i=1:n+1
Lix=1.0;
for j=1:n+1
if j ~= i
Lix=Lix.*(xint-x(j))/(x(i)-x(j));
end
end
yint=yint+y(i).*Lix;
end
end
function xp = GLL_nodes(N)
syms x
Lo = (1-x^2)*diff(legendreP(N,x),x);
xp = double(vpasolve(Lo == 0));
end