Sorry, I don't use MATLAB.
Your situation is quite explanable. For a single pure tone with a whole number of cycles in the frame, the expecatation is that two bins, a conjugate pair, will be non-zero and the rest of the bins will be zero. The relationship between the signal definition and those bin values using a 1/N normalized DFT is this:
$$ S_n = M cos( f \cdot \frac{ 2\pi }{ N } n + \phi ) \tag 4 $$
(For unnormalized, like Matlab, multiply by $N$)
$$ Z_k = \frac{ M }{ 2 } e^{ i \phi } \, \text{ when } \, k=f \tag {19} $$
Using the terminology in my article:
The conjugate bin at $N-k$ will have the conjugate value, that is, the phase has the opposite sign.
Now, for the rest of the bins, the values you are getting are meaningless and due to numerical imprecision. Technically, they should be zero. In polar coordinates the value at the origin has no phase angle, or it can be any angle if you prefer.
Now, if you have a pure tone that is not a whole number of cycles in the frame, the phase angle for each bin can be found without calculating the DFT itself by using equation (25).
If you have a frequency that is near a whole number of cycles, the relationship between the phase in the signal definition and the phase of the DFT bin is given in my answer here:
FFT Phase interpretation of input signal with non-integer number of cycles in FFT window