Speed up the code with legendre or Gegenbauer polynomials combined with the symbolic function?

Dear users

I have below code with Symbolic variables to evaluate a function. The a_n,b_n,c_n and d_n are numerical arrays which are calculated by a code which is not shown. And the "z" and "rho" are symbolic function defined in main code. It runs very slowly, someone told me I should vectorized it. After I have done vectorization (the below code is after the vectorization), it was slow as well. Then I realised that the time consuming part is the Gegenbauer polynomials or legendre polynomials. When the order of the Gegenbauer polynomials is small like 10, the speed is ok, but it is supper slow when n is 50. My questions are

(1) If I swich from Gegenbauer polynomials to legendre polynomials (as they can conver to each other), will it run faster when the n is 50 or even bigger?

(2)If not, how can I run the higher order (higher than 50) of Gegenbauer polynomials or legendre polynomials with symbolic function?

function [str_func_syms] = str_func_syms(c,z,rho,a_n,b_n,c_n,d_n,N)
%UNTITLED2 Summary of this function goes here
%   Detailed explanation goes here

zeta=atanh(2*c*z/(z.^2+rho.^2+c.^2));
% sin(eta)=(rho/z)*sinh(zeta);

% mu=cos(eta)=cosh(zeta)-(c/z)*sinh(zeta);
mu=cosh(zeta)-(c/z)*sinh(zeta);

n = 0:N;
poly_syms = (a_n .* cosh((n-1/2)*zeta)...
      + b_n .* sinh((n-1/2)*zeta)...
      + c_n .* cosh((n+3/2)*zeta)...
      + d_n .* sinh((n+3/2)*zeta))...
      .* gegenbauerC(n+1,-1/2,mu);

poly_syms_sum=sum(poly_syms);

str_func_syms=(cosh(zeta)-mu).^(-3/2)*poly_syms_sum;



end