Plotting phase delay from first principles giving different result from phasedelay() function. Why?
The problem is the definition of the phase. The command angle() computes the phase $\phi(\omega)$ according to

$$H(e^{j\omega})=|H(e^{j\omega})|e^{j\phi(\omega)}\tag{1}$$

where $H(e^{j\omega})$ is the complex frequency response. At frequencies $\omega$ where the frequency response has zeros, the phase jumps by $\pi$. This is shown in the two left figures of the plot below. You can't compute the phase delay from $\phi(\omega)$ because of the jumps. Note that in your computation, the phase delay is correct until the first phase jump due to the first zero of the magnitude response.

The phase delay needs to be computed from a continuous phase $\varphi(\omega)$, which is implicitly defined by

$$H(e^{j\omega})=A(\omega)e^{j\varphi(\omega)}\tag{2}$$

where $A(\omega)$ is a real-valued but bi-polar smooth function. This continuous phase is shown in the top right figure below, and it is obviously linear. From this continuous phase, the phase delay can be computed by

$$\tau_{ph}=-\frac{\varphi(\omega)}{\omega}\tag{3}$$

and for the given linear-phase FIR filter the phase delay is constant:

$$\tau_{ph}=\frac{N-1}{2}\tag{4}$$

where $N$ is the filter length.

Note that the continuous phase $\varphi(\omega)$ cannot be obtained from $\phi(\omega)$ using the Matlab function unwrap(), because the latter just removes jumps by (multiples of) $2\pi$ due to phase ambiguity. Here we have to deal with actual phase jumps due to zeros of $H(e^{j\omega})$.