Bosonic qubits are a promising route to building fault-tolerant quantum computers on a variety of physical platforms. Studying the performance of bosonic qubits under realistic gates and measurements is challenging with existing analytical and numerical tools. We present a novel formalism for simulating classes of states that can be represented as linear combinations of Gaussian functions in phase space. This formalism allows us to analyze and simulate a wide class of non-Gaussian states, transformations and measurements. We demonstrate how useful classes of bosonic qubits-Gottesman-Kitaev-Preskill (GKP), cat, and Fock states-can be simulated using this formalism, opening the door to investigating the behaviour of bosonic qubits under Gaussian channels and measurements, non-Gaussian transformations such as those achieved via gate teleportation, and important non-Gaussian measurements such as threshold and photon-number detection. Our formalism enables simulating these situations with levels of accuracy that are not feasible with existing methods. Finally, we use a method informed by our formalism to simulate circuits critical to the study of fault-tolerant quantum computing with bosonic qubits but beyond the reach of existing techniques. Specifically, we examine how finite-energy GKP states transform under realistic qubit phase gates; interface with a CV cluster state; and transform under non-Clifford T gate teleportation using magic states. We implement our simulation method as a part of the open-source Strawberry Fields Python library.