Accurate predictions for the spread and evolution of epidemics have significant societal and economic impacts. The temporal evolution of infected (or dead) persons has been described as an epidemic wave with an isolated peak and tails. Epidemic waves have been simulated and studied using the classical SIR model that describes the evolution of susceptible (S), infected (I), and recovered (R) individuals. To illustrate the fundamental dynamics of an epidemic wave, the dependence of solutions on parameters, and the dependence of predictability horizons on various types of solutions, Paxson and Shen (2022c, IJBC)  proposed a Korteweg-de Vries (KdV)-SIR equation and obtained analytical solutions. I will first introduce the KdV-SIR equation and discuss its relationship to the non-dissipative Lorenz model. I then illustrate how the KdV-SIR equation and its analytical solutions were applied for the estimates of predictability horizons using the COVID-19 data from the USA. The KdV-SIR model is mathematically identical to the non-dissipative Lorenz 1963 model and the KdV equation in a traveling-wave coordinate. As a result, the dynamics of an epidemic wave and its predictability can be understood by applying approaches used in nonlinear dynamics, and by comparing the aforementioned systems. For example, a typical solitary wave solution is a homoclinic orbit that connects a stable and an unstable manifold at the saddle point within the I - I’ space. The KdV-SIR equation additionally produces two other types of solutions, including oscillatory and unbounded solutions. The analysis of two critical points reveals the features of solutions near a turning point. Using analytical solutions and hypothetical observed data, a simple formula was derived for determining a predictability horizon (i.e., the timing for the peak of an epidemic wave). The above procedures were then applied to the real world COVID-19 data from the USA. To illustrate the effectiveness of the procedures in real data application, raw data were first processed using a Fast Fourier Transform (FFT), the empirical mode decomposition (EMD) method or the 28-days running mean methods. For analysis, smooth data were then used to generate ensemble predictions.