A numerical method for solving systems of nonlinear one-dimensional balance laws, based on multivariate sigmoidal functions approximation, is developed. Constructive approximation theorems are first established in both, uniform and $$L^p$$ Lp norms. A priori as well as a posteriori error estimates are derived for the numerical solutions, when various kinds of sigmoidal functions, such as unit step, logistic and Gompertz functions, are chosen. The residual of the numerical method is also estimated. Numerical examples are given to test the performance of the algorithm, a comparison with the Godunov method is made concerning accuracy and computational cost. Finally, the numerical stability of the method is analyzed.