Gauss-Laguerre quadrature rule

% matlab script to derive the 2-point Gauss-Laguerre quadrature rule% and use it on an example% inelegantly set up equations from method of undetermined coefficients% and solveclear all close allformat longeq1 = 'w1*1 + w2*1 = 1';eq2 = 'w1*x1 + w2*x2 = 1';eq3 = 'w1*x1^2 + w2*x2^2 = 2';eq4 = 'w1*x1^3 + w2*x2^3 = 6';[w1,w2,x1,x2] = solve(eq1,eq2,eq3,eq4)pause% note: there are two solutions% we pick the one where x1 < x2[x1,index] = min(double(x1))w1 = double(w1(index))x2 = double(x2(index))w2 = double(w2(index))% define the integrandf = @(x) exp(x).*log(1+exp(-x));% use Gauss-Laguerre quadrature to approximate itglq = w1*feval(f,x1) + w2*feval(f,x2)% now let's find the exact answer and see what the error is ...syms xI = double(int(log(1+exp(-x)),0,Inf))error = I - glqpause% or we can estimate the neglected tail% trying integral(f,0,Inf) or integral(f,0,realmax)% won't work!f = @(x) log(1+exp(-x));Q10 = integral(f,0,10)Q20 = integral(f,0,20)Q30 = integral(f,0,30)% from here we see that we have convergence to 8 decimal places% we can also perform a change of variable to make the interval finiteF = @(t) log(1+t)./t;integral(F,0,1)% check for problems at t=0ezplot(F,0,1)integral(F,realmin,1)