@Dan I cannot help you much with that in MATLAB. These days I mostly use R for anything beyond basic statistics. For your reference, here is a Q/A about something what I believe is exactly similar to your situation:
http://stats.stackexchange.com/questions/36351/how-to-run-a-two-way-anova-with-a-random-variable-followed-by-pairwise-compariso/36358#36358

@Patrizia your problem with NaN and in general, missing data is beyond the scope of this code. Generally this is handled separately through what we call imputation methods. They range from very basic (and unsatisfactory) approaches like removing incomplete cases to more advanced approaches like multiple imputations. However, the point is that there is no single, universally best approach for the missing data problem.

@J The main difference between repeated measures and the n-way ANOVA is that the former can model within subject variations while the latter can only model between subject variations. Essentially repeated measures ANOVA is a small subset of linear mixed models. In general you cannot model mixed models with simple, n-way ANOVA. In SPSS as well as the newest version of MATLAB linear mixed modeling is supported.

@Justin: I'm struggling with the comment editor. A previous comment was magically deleted and the other one does not make much sense now. In short, I verified the code and I think it is fine. Also, I run the same data through R and got identical results. You can see the R output here: http://imgur.com/FHJc4.png

As I had this problem and took me a while to figure it out, I thought it may be useful to someone else.
For anova_rm to be compatible with multcompare function of MATLAB, we need a stat output. I read somewhere (I printed out the webpage about a year ago but I can't find it anymore) that the post-hoc for repeated measures anova is the same as independent measures anova. You just need to replace MS(within) with MS(error) and df(within) with df(error) and then use Tuckey's equation to find the critical value. Thus, as I work with anova_rm as a one-way measure, I added the following lines at the end of anova_rm and added "stat" variable to the output of the function. I can now feed this stat to multcompare and the function sees stat as the output of anova1 and estimates the critical value using Tuckey (if that is your selected ctype) in the same way it does for anova1. As I don't use it for two-way measures, I haven't checked how to define stat for that purpose but I am sure it shouldn't be hard now.

@Dan I cannot help you much with that in MATLAB. These days I mostly use R for anything beyond basic statistics. For your reference, here is a Q/A about something what I believe is exactly similar to your situation:
http://stats.stackexchange.com/questions/36351/how-to-run-a-two-way-anova-with-a-random-variable-followed-by-pairwise-compariso/36358#36358

Well, I couldn't find that specific website but here is another one. Please let me know if you found any cite-able references for this purpose.
http://web.mnstate.edu/malonech/Psy232/Notes/Repeated%20Measures%20ANOVA%20GW14.htm

As I had this problem and took me a while to figure it out, I thought it may be useful to someone else.
For anova_rm to be compatible with multcompare function of MATLAB, we need a stat output. I read somewhere (I printed out the webpage about a year ago but I can't find it anymore) that the post-hoc for repeated measures anova is the same as independent measures anova. You just need to replace MS(within) with MS(error) and df(within) with df(error) and then use Tuckey's equation to find the critical value. Thus, as I work with anova_rm as a one-way measure, I added the following lines at the end of anova_rm and added "stat" variable to the output of the function. I can now feed this stat to multcompare and the function sees stat as the output of anova1 and estimates the critical value using Tuckey (if that is your selected ctype) in the same way it does for anova1. As I don't use it for two-way measures, I haven't checked how to define stat for that purpose but I am sure it shouldn't be hard now.
[n,k]=size(X{1});
stats.gnames=num2str([1:k]');
stats.n=n*ones(1,k);
stats.source='anova1';
stats.means=y_j';
stats.df=n*k-n-k+1;
stats.s=sqrt(msR);

@Dan I cannot help you much with that in MATLAB. These days I mostly use R for anything beyond basic statistics. For your reference, here is a Q/A about something what I believe is exactly similar to your situation:
http://stats.stackexchange.com/questions/36351/how-to-run-a-two-way-anova-with-a-random-variable-followed-by-pairwise-compariso/36358#36358

Comment only