Hi,
The NY and NU sizes are based on the number IOs (inputs and outputs) your system has and would be fixed. You can compute the response at various frequencies with the resulting mag and phase being a NY-by-Nu-by-Nw vector where Nw is the number of frequency points. If you are trying to extract data at more frequencies or a specific set of frequencies you can specify the frequency vector explicitly
[mag,phase,wout] = bode(sys,win)
In this case win = wout. For example, evaluating 100 points equally logrithmic spaced points between 10^-1 and 10 would be done has follows
[mag,phase,wout] = bode(sys, logspace(-1,1,100))