I want to extract w and x from this code, but unfortunately i cant and it gives me only 1 answer.

function [sol]=numintla(fun,n)
[xx,w]=gaussla(6); % Call the gaussla.m file
sx=size(xx,1);
sol = 0;
for i=1:sx,
    x=xx(i);
    fx=eval(fun);
    sol=sol+w(i)*fx;
end
end
function [xx,ww]=gaussla(n)
format long;
b(1)=1;
for i=1:n, % Calculation of the binomial coefficients
    b(i+1)=(factorial(n))/(factorial(i)*(factorial(n+1-i)));
end
for i=1:n+1, % The polynomial coefficients
    poly(i)=((-1)^(n+1-i))*b(i)*(factorial(n)/(factorial(n+1-i)));
end
xx=roots(poly);% The polynomial roots
for i=1:n, % Coefficients of the first derivative of the polynomial
    polycd(i)=poly(i)*(n+1-i);
end
for i=1:n, % Evaluation
    x=xx(i);
    solde=0;
    for k=1:n,
        solde=solde+polycd(k)*(x^(n-k));
    end
    ww(i,1)=(factorial(n)^2)/(xx(i)*(solde^2));
end
end

Answer 1

[Walter Roberson]

function [sol, w, x] = numintla(fun,n)

I would recommend that you reconsider your requirement. Your x variable is whatever was last written into x by your for i loop. It would make more sense if you were to return xx instead of x.

Answer 2

[Riccardo Scorretti]

Hi Alireza,

I'm afraid that the fact your code returns only one value is the last of your problems. More specifically, you provides two functions: one returns two values and the other returns two. So, it depends on what you want to do: if your purpose is to get the quadrature formula, you can invoke guass1a in this way (the number 3 is just an example: you can use any positive integer):

[x, w] = gaussla(3)

In my opinion (but I may be wrong) the true problems are more fundamental. You seems to be willing to compute (gaussla) the nodes and the weights of a Gauss-Legendre quadrature formula, which is usually defined in the range [-1 1]. The function numintla integrates numerically the function fun over a not precised domain (which I suppose is [-1 1]).

I'm afraid the function gaussla returns the wrong result. For instance, you can check the first coefficients here: https://en.wikipedia.org/wiki/Gaussian_quadrature. For n = 1 and n = 2 you get respectively:

x =

1

w =

1

and

x =

1.000000000000000 + 1.000000000000000i

1.000000000000000 - 1.000000000000000i

w =

-0.500000000000000 + 0.500000000000001i

-0.500000000000000 - 0.500000000000001i

One observes that coefficients are complex (this is wrong). The error is likely ot be in the computation of the coefficients of Legendre polynomials (see https://en.wikipedia.org/wiki/Legendre_polynomials).

If I'm not wrong, I suggest you to have a look this code: https://fr.mathworks.com/matlabcentral/fileexchange/4540-legendre-gauss-quadrature-weights-and-nodes