Based on the Green's functions that reflect mathematical properties of partial differential equations (PDE), we developed a novel preconditioner using neural networks (NNs) with high accuracy and small computational cost for improving the convergence property of an iterative implicit solver. As the dense and uniform computation involved in NNs are more efficient than that of the conventional PDE solver schemes, we could solve the time evolution of a 405,017,091 degrees-of-freedom highly heterogeneous problem in 5.48-fold shorter time compared to a typical PDE solver. The method is also suitable for use with low-precision arithmetic in NNs as the accuracy of the final solution is guaranteed. The localized property of NNs enable high scalability for solving large problems (103,305,758,211 degrees-of-freedom problem solved with 97.4% weak scalability using 256 Cascade Lake Xeon CPU-based Oakbridge-CX nodes with a total of 14336 CPU cores with developed MPI-OpenMP hybrid code). This method can be used in various PDE-based simulations and has potential to make broad ripple effects in various fields.