syms zeta mu rho z eta
zeta=atanh(2*c.*z./(z.^2+rho.^2+c.^2));
cos(eta)=cosh(zeta)-(c/z)*sinh(zeta);
mu=cos(eta);
for n=0:N
poly_syms(n+1,:) =(a_n(n+1,:)*cosh((n-1/2)*zeta)...
+b_n(n+1,:)*sinh((n-1/2)*zeta)...
+c_n(n+1,:)*cosh((n+3/2)*zeta)...
+d_n(n+1,:)*sinh((n+3/2)*zeta))...
.*gegenbauerC(n+1,-1/2,mu);
end
poly_syms_sum=sum(poly_syms);
str_func_syms=(cosh(zeta)-mu).^(-3/2)*poly_syms_sum;
poly_func=matlabFunction(poly_syms);
str_func_func=matlabFunction(str_func_syms);
[domain_rho,domain_z]=meshgrid(rho_interval,z_interval);
str_func=str_func_func(domain_rho,domain_z);
