In this paper, we describe a numerical approach based on Bernstein matrix method to solve a mathematical model arising from a model of neuronal variability. For this purpose, first we present some preliminaries about Bernstein polynomials then a new approach implementing a collocation method in combination with matrices of Bernstein polynomials is introduced to approximate the solution of these equations with variable coefficients under the boundary conditions. Finally, a test example is solved to demonstrate the efficiency of the presented numerical scheme.