A numerical model of Gaussian beam propagation in planar Fabry-Perot (FP) etalons is presented. The model is based on the ABCD transfer matrix method. This method is easy to use and interpret, and readily connects models of lenses, mirrors, fibres and other optics to aid simulating complex multi-component etalon systems. To validate the etalon model, its predictions were verified using a previously validated model based on Fourier optics. To demonstrate its utility, three different etalon systems were simulated. The results suggest the model is valid and versatile and could aid in designing and understanding a range of systems containing planar FP etalons. The method could be extended to model higher order beams, other FP type devices such as plano-concave resonators, and more complex etalon systems such as those involving tilted components.