The phase velocity in your case should not vary much, so i guess that fitting the CFL to the theoreticaly fastest wave should work fine. If you need an accurate and least dispersive method than implicit methods are not very good since they tend to damp the solution so i'd go with an explicit scheme. If you choose FD method then you should be aware that Arakawa A, B, C, D and E grids (staggering methods) produce dispersion (even for the waves that are nondispersive), also you have to be aware that choosing between leepfrog, forward and backward in time produces different dispersions and can affect stability as well. So, what you need is to idealize the problem and solve it analiticaly, and then to employ a couple of methods on this idealized problem, and see which one works best, and then to use this method on nonidealized - real problem, and cross your fingers and hope it'll work. This is really the way to do it.
Coming back to the varying phase velocity, this should not be a problem, in meteorology this is the problem that is being solved every day since atmospheric motion is really a wave motion (atmospheric equations of motion can be transformed with some aproximations to the Helmholtz eq), and there are many phase speeds that are involved - from ~ 1cm/s up to the speed of sound. There is no numerical method that does not distort any wave, every method in the book _WILL_ distort your waves, it is up to you to find the method whose distortion will be most acceptable for your problem.
The problem you encounter is well known problem in meteorology. The two methods that are being used most are FE and spectral method. The former is simpler and more flexible and the later being more accurate since it keeps the derivations analytic. The problem with the spectral method is that you have to make your problem periodic, and you have to go from grid-point space to spectral space every timestep. In regard to stability, this is not a problem in many FE methods as long as you fulfill Courant-Friedrich-Levi criterion i.e. as long as your timestep is shorter than the time needed for the fastes wave to pass the distance between two grid points. When aplying FE method you have to make a choice of a scheme and depending on the scheme you can expect different distortions of your wave.
Try to look in the book Numerical Methods for Solving Wave Equiations in Geophysical Fluid Dynamics by Dale Durran. Good on-line source is http://www.ecmwf.int/newsevents/training/lecture_n otes/LN_NM.html especially "Numerical methods" by R. W. Riddaway and M. Hortal
The phase velocity in your case should not vary much, so i guess that fitting the CFL to the theoreticaly fastest wave should work fine. If you need an accurate and least dispersive method than implicit methods are not very good since they tend to damp the solution so i'd go with an explicit scheme. If you choose FD method then you should be aware that Arakawa A, B, C, D and E grids (staggering methods) produce dispersion (even for the waves that are nondispersive), also you have to be aware that choosing between leepfrog, forward and backward in time produces different dispersions and can affect stability as well. So, what you need is to idealize the problem and solve it analiticaly, and then to employ a couple of methods on this idealized problem, and see which one works best, and then to use this method on nonidealized - real problem, and cross your fingers and hope it'll work. This is really the way to do it. Coming back to the varying phase velocity, this should not be a problem, in meteorology this is the problem that is being solved every day since atmospheric motion is really a wave motion (atmospheric equations of motion can be transformed with some aproximations to the Helmholtz eq), and there are many phase speeds that are involved - from ~ 1cm/s up to the speed of sound. There is no numerical method that does not distort any wave, every method in the book _WILL_ distort your waves, it is up to you to find the method whose distortion will be most acceptable for your problem.
The problem you encounter is well known problem in meteorology. The two methods that are being used most are FE and spectral method. The former is simpler and more flexible and the later being more accurate since it keeps the derivations analytic. The problem with the spectral method is that you have to make your problem periodic, and you have to go from grid-point space to spectral space every timestep. In regard to stability, this is not a problem in many FE methods as long as you fulfill Courant-Friedrich-Levi criterion i.e. as long as your timestep is shorter than the time needed for the fastes wave to pass the distance between two grid points. When aplying FE method you have to make a choice of a scheme and depending on the scheme you can expect different distortions of your wave. Try to look in the book Numerical Methods for Solving Wave Equiations in Geophysical Fluid Dynamics by Dale Durran. Good on-line source is http://www.ecmwf.int/newsevents/training/lecture_n otes/LN_NM.html especially "Numerical methods" by R. W. Riddaway and M. Hortal