In the absence of inhibition, excitatory neuronal networks can alternate between bursts and interburst intervals (IBI), with heterogeneous length distributions. As this dynamic remains unclear, we develop here a bursting model based on synaptic depression and facilitation that also accounts for afterhyperpolarization (AHP), which is a key component of IBI. The framework is a novel stochastic three dimensional dynamical system perturbed by noise: numerical simulations can reproduce a succession of bursts and interbursts exhibiting AHP refractory periods [1] . The Matlab source code is organized in three sub folders: 01_Functions Contains the code to simulate the model, the function used to segment the time series into Burst, AHP and QP phases, the code used to obtain the characteristics and the analytical formulas for burst and AHP durations presented in [1]. 02_Figures Contains one script for each figure in [1], these script can be run directly from the root of the code folder and will load the correct parameter values and data (if necessary) to reproduce the figures presented in the article. 03_Data Contains the simulated data and parameter values. Please quote the zenodo or doi reference in publications when using the code.