function x = legendreRoots(n,method)

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,:));
    if strcmp(method,'newton')
        x = x - fx./dfx;
    elseif strcmp(method,'halley')
        d2fx = (n-1)*n^2*(n+1)/((2*n+1)*(2*n-1))*1./(1-x.^2).^2 .* (Px(n-1,:) - Px(n+1,:)) ...
            - n*(n+1)^2*(n+2)/((2*n+1)*(2*n+3))*1./(1-x.^2).^2 .* (Px(n+1,:) - Px(n+3,:)) ...
            + n*(n+1)/(2*n+1) * 2*x./(1-x.^2).^2 .* (Px(n,:) - Px(n+2,:));
        x = x - 2*fx.*dfx./(2*dfx.^2-fx.*d2fx);
    else
        fprintf(2,'Error: Method is unkown. Use method = {newton,halley}.\n')
        return
    end   
    
    %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