Exact solutions for one-dimensional steady-state and transient liquid flow toward a water table in heterogeneous unsaturated soils are critical in predicting saturation profiles in several real-world applications including interpretation of climate change effects on the subsurface and impacts on slope stability. In this study, vertical heterogeneity in saturated hydraulic conductivity with depth is characterized by an exponential decay function. A steady-state solution is derived based on Darcy's law and the water table depth, and two transient-state solutions are obtained from Richards' equation using the Laplace transform and the modified Bessel equation under common upper boundary conditions, that is, flow rate and pressure head, following the initial steady-state flow condition and a water table at a specified depth. The transient and steady-state solutions are compared with numerical solutions obtained from a multi-layered configuration and an analytical solution for homogeneous soils, demonstrating their reliability and efficacy. Various hypothetical heterogeneous soils with differing parameters are employed to assess the flow behaviors, illustrating that vertical heterogeneity impacts the pressure head profiles. More pronounced heterogeneous soils exhibit lower pressure head and effective saturation during flow compared to homogeneous soils, which depend on the air-entry value and the water table depth, irrespective of the upper boundary conditions. The solution has been employed to predict volumetric water content profiles in the Loess Plateau of China and evaluate the effect of infiltration on the stability of a shallow landslide in Japan based on the infinite slope model and the suction stress concept.