We have examined the validity of the time-dependent variational approximation (TDVA) to the Gaussian wavepacket method (GWM) for quantum double-well (DW) systems, by using the quasiexact spectral method (SM). Comparisons between results of wavefunctions, averages of position and momentum, the auto-correlation function, and an uncertainty product calculated by SM and TDVA have been made. It has been shown that a given initial Gaussian wavepacket in SM is quickly deformed at t > 0 where a wavepacket cannot be expressed by a single Gaussian, and that assumptions on averages of higher-order fluctuations in TDVA are not justified. These results cast some doubt on an application of TDVA to DW systems. Gaussian wavepacket dynamics in anharmonic potential systems is studied also.