The proper orthogonal decomposition (POD) method for analysis of nonlinear panel flutter subjected to supersonic flow is presented. Optimal POD modes are extracted from a chaotic Galerkin mode responses. The aeroelastic equations of motion are constructed using von Karman plate theory, first-order piston theory and quasi-steady thermal stress theory. A simply-supported plate with thermal loads from a uniformly distributed temperature is considered. Many types of panel behaviors, including stable flat, dynamically stable buckled, limit cycle oscillation, nonharmonic periodic motion, quasi-periodic motion and chaotic motion are observed. Our primary focus is on chaos and the route to chaos. It is found that a sudden transition from the buckled state to chaos occurs. Time history, phase portrait, Poincaré map, bifurcation diagram and Lyapunov exponent are employed to study chaos. The POD chaotic results obtained are compared with the traditional Galerkin solutions. It is shown that the POD method can obtain accurate chaotic solutions, using fewer modes and less computational effort than the Galerkin mode approach; additionally, the POD method converges faster in the analysis of chaotic transients. Effects of length-to-width ratios and thermal loads are presented. It is found that a smaller width for fixed length will produce more stable flutter response, while the thermal loads degrade the flutter boundary and result in a more complex evolution of dynamic motions. The numerical simulations show that the robustness of the POD modes depends on the dynamic pressure but not on temperature.