The puncture method for dealing with black holes in the numerical simulation of vacuum spacetimes is remarkably successful when combined with the BSSN formulation of the Einstein equations. We examine a generalized class of formulations modeled along the lines of the Laguna-Shoemaker system, including BSSN as a special case. The formulation is a two parameter generalization of the choice of variables used in standard BSSN evolutions. Numerical stability of the standard finite difference methods is proven for the formulation in the linear regime around flat space, a special case of which is the numerical stability of BSSN. Numerical evolutions are presented and compared with a standard BSSN implementation. We find that a significant portion of the parameter space leads to stable evolutions and that standard BSSN is located near the edge of the stability region. Non-standard parameter choices typically result in smoother behaviour of the evolution variables close to the puncture and thus hold promise for improved accuracy in, e.g., long-term BH binary inspirals, and for overcoming (numerical) stability problems still encountered in some types of black-hole simulations, e.g., in $D \ge 6$ dimensions.