chebyshev spline and pchip usage

Dear Matlab Users, I want to use equally spaced nodes or chebyshev nodes to evaluate my function at the nodes and then I want to construct an interpolating polynomial. Then I want to find zeros of the interpolating polynomial in an interval Z=[-1,1]

kd=constant; lambda=constant; Delta_E=constant;

z_dot= -((((1+z)./2)).*exp(1).^-(kd*(sqrt(2*lambda*(1+z)))))*(1+kappa)*(sqrt(1-z^2))*(sin(phi))/2;

Phi_dot=((Delta_E)+(lambda*z)-(((sqrt((1+z)./2)).*exp(1).^(-(kd*(sqrt(2*lambda*(1+z))))))).*(kappa*(sqrt((1-z)./(1+z)))-(sqrt((1+z)./(1-z)))))*(cos(phi))/2;

I would highly appreciate if you let me know how I can start, then I will write the code and ask again if I stuck somewhere.

Best, Fatemeh