Density estimation is hard.
With the CDF, we can always use the Empirical Distribution Function (EDF) to get a consistent estimate of the CDF under mild assumptions (random independent draws from target distribution).
This is definitively not the case for density estimation. There is no way to get good estimates of the PDF without either enforcing smoothness (kernel density estimation or log-splines, for example) or other assumptions (belongs to parametric family, is log-concave, etc). This is because without any smoothness assumptions/priors, the distribution that would make the data we have seen most likely would be the distribution that puts probability $1/n$ on each of the data points we observed. Since we observed a finite data set, this is necessarily a discrete random variable, for which the PDF is undefined (i.e., it is defined by a probability mass function, not probability density function).
So let's dive into smoothness and assumptions. Parametric assumptions will generally have the best operating characteristics (i.e. lowest MSE) when the assumptions hold true, but can be heavily biased when this is false. As such, these methods are usually discarded for density estimation; choosing the wrong parametric family will have serious consequences (unlike, say linear regression, in which case OLS is known to still be a reasonable choice even if the data is not normal), and we typically have no good reason to believe that we know exactly what family the distribution belongs to.
If we decide to reject parametric assumptions, we can move to non-parametric models. Three of the most popular methods are kernel density estimation, Gaussian mixture models and log-splines. I do believe kernel density estimation is probably considered most popular, although that may differ by application. In all of these, we will need to select some sort smoothing parameter (bandwidth for kernel density estimation, number of components for GMM and number of knots for log-splines). Selection of smoothing parameter is crucial, yet non-trivial. There are automatic methods for selection (e.g., AIC), but be warned that these methods can still lead to some unsatisfying results.
In my experience, I've found kernel density estimation to be most reliable. In fact,
density function uses this and the default settings tend to be quite reasonable. But be warned that you don't have to look too hard to find cases where the results may be unsatisfactory.