Chebyshev Economization on Fresnel Cosine Function.

I wrote the following function to calculate the Chebyshev Economization of the Fresnel Cosinus Function.

The first while loop will calculate my taylor series and the second while loop will keep adding chebyshev polynomials until the approximation is accurate enough.

if true
  function [] = taylorFresnelCFunction()
  xx = 0:0.05:5;
  for i =1:length(xx)
      y(i) = fresnel(xx(i));
  end
  n = 0;
  y_taylor = @(x) (-1)^n*(pi/2)^(2*n)*x.^(4*n+1)/(factorial(2*n)*(4*n+1));
  yApprox = y_taylor(xx);
  difference = max(abs(y-yApprox));
  while difference >= 0.05
      n = n + 1;
      temp = @(x) (-1)^n*(pi/2)^(2*n)*x.^(4*n+1)/(factorial(2*n)*(4*n+1));
      tEF = @(x) ( y_taylor(x) + temp(x) );
      difference = max(abs(tEF(xx)-y));
      y_taylor = tEF;
  end
  f = tEF(xx);

    differenceChebyshev = 1;
    n = 1;
    while differenceChebyshev >= 0.01
        T = zeros(n,n);
        x = zeros(1,n);     % Maak het polynoom f(x) = x aan.
        x(2) = 1;
        T(1,1) = 1;         % T_0 = 1
        T(2,2) = 1;         % T_1 = x
        for i = 3:n
            temp = conv(x,T(i-1,:));
            T(i,:) = 2*temp(1:n) - T(i-2,:);    % T_n = 2*x*T_(n-1)-T_(n-2)
        end
        T = fliplr(T);

        D = zeros(length(xx),n);

        for i=1:length(xx)
            for j=1:n
                D(i,j) = polyval(T(j,:),xx(i));
            end
        end
        coefficients = D\f'
        chebyshevPolynomial = bsxfun(@times,coefficients,T);
        chebyshevPolynomial = sum(chebyshevPolynomial,1);
        chebPolDer = polyder(chebyshevPolynomial);

        x1 = 0:0.05:5;
        y1 = polyval(chebyshevPolynomial,x1);

        differenceChebyshev = max(abs(y1-f));
        n = n + 1;
    end
    n = n - 1;
   hold all
   plot(xx,y1)
   plot(xx,y)
end
end

My problem is that I already wrote this function for the Error Function (taylorErrorFunction) and this works without a problem.

I thought replacing the taylor series in taylorErrorFunction with the appropriate series for the Fresnel Cosine would give me the result but this didn't happen. I'm getting 'NaN' as result for the coefficients.

Does anyone know why the fresnel cosine is not being calculated correctly? I have added the taylorErrorFunction code (so you can see that it does work for that particular function) as well as the file written to calculate fresnel cosine (and sine) in case anyone wants to run the program.

I am using the Fresnel Cosine with the argument (pi/2*t²).