An efficient second-order stable numerical method
is presented to solve the model partial differential equations
of thermal glass fiber processing.
The physical process and structure of the model equations are described first.
The numerical issues are then clarified.
The heart of our method is a MacCormack scheme with flux limiting.
The numerical method is validated on a linearized isothermal model
and by comparison with known exact stationary solutions.
The numerical method is then generalized to solve the equations
of motion of thermal glass fiber drawing, exhibiting order of convergence.
Further, the nonlinear PDE scheme is benchmarked against
an independent linearized stability analysis of boundary value solutions
near the onset of instability, which demonstrates the efficiency of the method.