Regarding Shuo Han's alternative solution. For robustness I had to modify it to the below. My theory is that the division with V sometimes introduces numerical errors that are large enough to result in a new negative eigenvalue.

[V,D] = eig(A);
A_psd = V * diag(max(diag(D),eps)) / V;

Regarding Shuo Han's alternative solution. For robustness I had to modify it to the below. My theory is that the division with V sometimes introduces numerical errors that are large enough to result in a new negative eigenvalue.
[V,D] = eig(A);
A_psd = V * diag(max(diag(D),eps)) / V;

Comment only