We present a numerical algorithm to implement entropy-based ($M_N$) moment closures in the context of a simple, linear kinetic equation that models particles moving through a material slab. The algorithm has two main, coupled components: a second-order kinetic scheme to update the PDE and a Newton-based solver for the dual of the optimization problem that defines the closure.
We study in detail the difficulties of solving the dual problem near the boundary of realizable moments, where quadrature formulas are less reliable and the Hessian of the dual objection function is highly ill-conditioned. Extensive numerical experiments are performed to illustrate these difficulties. In cases where the dual problem becomes ``too difficult'' to solve numerically, we propose a regularization technique to artificially move moments away from the realizable boundary in a way that still preserves local particle concentrations. We present results of numerical simulations for two challenging test problems in order to quantify the characteristics of the optimization solver and to investigate when and how frequently the regularization is needed.