I find the same thing. I got 52/1000 rejected the first way. That's close to the 5% value that you would expect. I got fewer the second way. That's because by estimating the parameter from the data used in the test, you are using a fit that is "too good." Imagine if you estimated the p parameter of a binomial distribution -- this would always give an exact match between the observed and expected counts of 1's and 0's. This case is not so extreme, but is similar.
The chi2gof function does provide a way to do an approximate adjustment for that, if you tell it that you have estimated one parameter:
[H(sim) P(sim) STATS] = chi2gof(X,'cdf',@(z)poisscdf(z,mean(X)),'nparams',1);
The concept is described further in lecture 12 of the MIT course notes.
As for your other question, the chi2gof function defines a set of bin edges and calls your function to compute the cdf values at these edges. The "@(z)" indicates that your function accepts one input to be represented by the symbol z, and the rest of the expression defines how the cdf is to be computed at the provided z values.
