function [xi,w]=GaussByNewtonMaehly(a,b,n);
% GAUSSBYNEWTONMAEHLY Gauss nodes and weights using Newton Maehly
%   [xi,w]=GaussByNewtonMaehly(a,b,n) computes the n nodes and weights
%   for Gauss Legendre quadrature for the interval (a,b). The nodes
%   are computed with the Newton-Maehly algorithm, and the weights by
%   using the function Phi and the three term recurrence relation.
%   Since the rule is symmetric, only half of the nodes and weights
%   are computed.

xi=[];
for k= 1:n                          % coefficients for Gauss-Legendre 
  beta(k)=k/sqrt(4*k^2-1);          % alpha(i)=0
end
phi=sqrt(3/2)*(b-a);                % =p_1'*int_a^b
p=1/sqrt(2);                        % =p_0
xm=(a+b)/2; xl=(b-a)/2;             % change of variables
anz=fix((n+1)/2); x0=1; 
for k=1:anz,
  y=x0; finished=0;
  m=2;                              % Newton double step
  while ~finished,
     x=y;                           % evaluate polynom and derivative
     p1s=0; p1=p; p2s=sqrt(3/2); p2=sqrt(3/2)*x;
     for i=2:n
       p0s=p1s; p0=p1; p1s=p2s; p1=p2;
       p2s=(p1+x*p1s-beta(i-1)*p0s)/beta(i);
       p2=(x*p1-beta(i-1)*p0)/beta(i);
     end
     s=sum(1./(x-xi(1:k-1))+1./(x+xi(1:k-1)));
     y=x-m*p2/p2s/(1-p2/p2s*s);     % Newton-Maehly step
     if y>=x,
        if m==1, finished=1;
        end
        if ~finished;               
          x0=x; m=1;                % stop double step
          y=x-p2/p2s/(1-p2/p2s*s);  % Neton backstep
        end
     end
  end
  xi(k)=y;
  phi1=0; phi2=phi; 
  for i=2:n
    phi0=phi1; phi1=phi2; phi2=(y*phi1-beta(i-1)*phi0)/beta(i);
  end
  w(k)=xl*phi2/p2s; w(n+1-k)=w(k);
end
for k=1:anz,                        % backchange variables
  y=xi(k); xi(k)=xm+xl*y; xi(n+1-k)=xm-xl*y;
end
xi=xi(:); w=w(:);