I am not sure about what is happening. Probably a problem with the computation of the binomial coefficients. These must be computed recursively, because the gamma computation is bad for high values of the argument. Remember that the binomial coefficients are equal to (-d)k/k! where I represented the Pochammer coefficient by (-d)k. Besides, if you are using Matlab, I suggest you to use other variables, because i and j are used in Matlab for representing the imaginary unity sqrt(-1). Try this. I would you suggest that, at the end, remove the initial calculated points because they are bad due to the transient of the derivative operator.