Partial differential equations (PDEs) fitting scientific data can represent physical laws with explainable mechanisms for various mathematically-oriented subjects. Most natural dynamics are expressed by PDEs with varying coefficients (PDEs-VC), which highlights the importance of PDE discovery. Previous algorithms can discover some simple instances of PDEs-VC but fail in the discovery of PDEs with coefficients of higher complexity, as a result of coefficient estimation inaccuracy. In this paper, we propose KO-PDE, a kernel optimized regression method that incorporates the kernel density estimation of adjacent coefficients to reduce the coefficient estimation error. KO-PDE can discover PDEs-VC on which previous baselines fail and is more robust against inevitable noise in data. In experiments, the PDEs-VC of seven challenging spatiotemporal scientific datasets in fluid dynamics are all discovered by KO-PDE, while the three baselines render false results in most cases. With state-of-the-art performance, KO-PDE sheds light on the automatic description of natural phenomenons using discovered PDEs in the real world.