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 behavior 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 continuous-variable 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.