I notice two problems in your code:
- at $f=500\,\text{kHz}$ you require a phase delay of $107$ samples, but your filter length is only $50$ taps; you can't have a delay that's greater than the filter length of an FIR filter.
- it looks like you normalize all frequencies by the sampling frequency; however, in Matlab a normalized frequency $f=1$ corresponds to the Nyquist frequency, i.e., $f_s/2$, not to $f_s$.