function x = rootsLegendre(n)

tol = 1e-12;
maxit = 1000;

% Check trivial cases
if n == 0
    x = [];
    return
elseif n == 1
    x = 0;
    return
elseif n == 2
    x(1) = -1/sqrt(3);
    x(2) = 1/sqrt(3);
    return
end

% Set initial guess
x = ones(1,n)*(1-(n-1)/(8*n^3)).*cos((4*(1:n)-1)/(4*n+2)*pi);
m = length(x);

%set inital value and calculate Pn(x0)
Px = zeros(n+3,m);
Px(1,:) = ones(1,m);
Px(2,:) = x;
for k = 1:n+1
    Px(k+2,:) = 1/(k+1)*((2*k+1)*x.*Px(k+1,:)-k*Px(k,:));
end
Pn = Px(n+1,:);

%*** perform fixpoint iteration
nit = 1;
while max(abs(Pn)) > tol && (nit <= maxit)
    fx = Pn;
    dfx = n*(n+1)/(2*n+1)*1./(1-x.^2).*(Px(n,:)-Px(n+2,:));
    % Newton step 
    x = x - fx./dfx;
    
    % calculate Pn(x)
    Px(2,:) = x;
    for k = 1:n+1
        Px(k+2,:) = 1/(k+1)*((2*k+1)*x.*Px(k+1,:)-k*Px(k,:));
    end
    Pn = Px(n+1,:);
    nit = nit + 1;
end

%*** transform to [-3,3]
x = 3*x;
