Error in power spectral density computation
Matlab splits your signal into 8 overlapping segments (by default). The difference I see between Matlab and your algorithm is that Matlab also multiplies each segment by a Hamming window. You should get pretty close results if you do that.