I believe you just need to redefine the critfun function from the one in the example:
function dev = critfun(X,Y)
model = fitglm(X,Y,'Distribution','binomial');
dev = model.Deviance;
replacing the critical value with
You might want to rename that variable something like rsqr, to avoid confusion.
After reading that example, and thinking about it a bit more, there might be some other nuances. That example states, "Adding a feature with no effect reduces the deviance by an amount that has a chi-square distribution with one degree of freedom". I'm not sure the same is true for R^2. So, that might bear some thought.
Also, I believe the deviance measure is something that is minimized, whereas R^2 is maximized. There is probably an adjustment that needs to be made for that as well. (One simple possibility would be to return 1-R^2 in the critical function, I guess.)