Plotting the zeros of jacobi polynomials

I would like to plot the zeros of jacobi polynomials of increasing degree, for which I have made the two functions at the bottom of this question.

Whenever I try to run "jacobiroots(50,1,2,3)", I get the error "Undefined function 'sym2poly' for input arguments of type 'function_handle'. Error in jacobiroots (line 7) coeff=sym2poly(p);". It's been a long time since I've worked with MATLAB so I don't really see how to solve this.

Thank you in advance!

function jacobiroots(n,a,b,c)
    syms x
    for k=1:n
        p=jacobirecursive(a*k,-(a+b)*k,(b+c)*k);
        coeff=sym2poly(p);
        r=roots(coeff);
        scatter(real(r),imag(r),'filled','red');
        title(k);
        pause(0.3)
    end
end

function output = jacobirecursive(n,alfa,beta)
    if n==1
        output = @(x) (alfa+1)+(alfa+beta+2)*(x-1)/2;
    else 
        first = (2*n+alfa+beta-1)*((2*n+alfa+beta)*(2*n+alfa+beta-2)*x+alfa^2-beta^2);
        second = 2*(n+alfa-1)*(n+beta-1)*(2*n+alfa+beta);
        denom = 2*n*(n+alfa+beta)*(2*n+alfa+beta-2);
        output = @(x) (first*jacobirecursive(n-1,alfa,beta)-second*jacobirecursive(n-2,alfa,beta))/denom;
    end
end

Answer 1

[Torsten]

I'd use

https://de.mathworks.com/help/symbolic/jacobip.html

:-)

[Catinca Mujdei]

I used this function before, but it does not give the desired accuracy; the zeros of the jacobi polynomials will tend to move away from the curve that they should be converging to. That is why I implemented the jacobi polynomials manually with a resursive function, hoping to afterwards be able to change the precision with the digits command.

[Torsten]

I don't have MATLAB available, but this might work:

function jacobiroots(n,a,b,c)
    for k=1:n
        p=jacobirecursive(a*k,-(a+b)*k,(b+c)*k);
        coeff=sym2poly(p);
        r=roots(coeff);
        scatter(real(r),imag(r),'filled','red');
        title(k);
        pause(0.3)
    end
end

function output = jacobirecursive(n,alfa,beta)
    syms x
    if n==1
        output = (alfa+1)+(alfa+beta+2)*(x-1)/2;
    else 
        first = (2*n+alfa+beta-1)*((2*n+alfa+beta)*(2*n+alfa+beta-2)*x+alfa^2-beta^2);
        second = 2*(n+alfa-1)*(n+beta-1)*(2*n+alfa+beta);
        denom = 2*n*(n+alfa+beta)*(2*n+alfa+beta-2);
        output = (first*jacobirecursive(n-1,alfa,beta)-second*jacobirecursive(n-2,alfa,beta))/denom;
    end
end

[Catinca Mujdei]

Thank you! It works. I also forgot to add an if statement to the function jacobirecursive in case n=0, which I now solved.