I think there's a bug in the code. if I see correctly the sequential nature of the testing will only be reflected if you set the adjusted p-value to the max across the current adjusted p and *all* preceding ones (as e.g. indicated in the respective wiki entry - the original paper has nothing on adjusted p-values afais).
now,in line 112 your algorithm picks the max of a given p-value and its immediate predecessor. this would be fine if it was done in a sequential fashion (each column updated in turn). but the way it stands the results can differ from adjusting to the max between all p-values up to the current one. you can end up with a vector of adjusted p-values that has decreasing p-values and significant results although the first test was non-significant (try e.g. [p,h]=bonf_holm([.01:.001:.02])). so in the current version non-significant values will be flagged as significant. a remedy would be to use a loop (ugly, I know) like:
also, you might want to add a
to avoid p-values >1