How to compute correct phase in FFT even after applying phase unwrap and zeroing round off error?
There's an error in your script. Your sampling frequency is not 100 Hz like you intended but is actually 99.9 Hz.  Therefore you don't actually have a whole number of cycles in your FFT which creates some leakage. This leakage will also impact the phase estimation.
If you want 1000 points at 100 Hz, the last point should not end at 10 seconds but at 10 seconds - 1/100 = 9.99. You should be careful when using linspace to define a time vector.
Use this n = 1 : 1000; t = (n-1)/100;