In the present study, three dimensional (3D) numerical modeling strategies of a thermosetting pultrusion process are investigated considering both transient and steady state approaches. For the transient solution, an unconditionally stable alternating direction implicit Douglas-Gunn (ADI-DG) scheme is implemented as a first contribution of its kind in this specific field of application. The corresponding results are compared with the results obtained from the transient fully implicit scheme, the straightforward extension of the 2D ADI and the steady state approach. The implementation of the proposed approach is described in detail. The calculated temperature and cure degree profiles at steady state are found to agree well with results obtained from similar analyses in the literature. Detailed case studies are carried out investigating the computational accuracy and the efficiency of the 3D ADI-DG solver. It is found that the steady state approach is much faster than the transient approach in terms of the computational time and the number of iteration loops to obtain converged results for reaching the steady state. Hence, it is highly suitable for automatic process optimization which often involves many design evaluations. On the other hand sometimes the transient regime may be of interest and here the proposed ADI-DG method shows to be considerably faster than the transient fully implicit method which is generally used by the general purpose commercial finite element solvers. Finally, using the proposed steady-state approach, a design of experiments is carried out for the curing characteristic of the product based on pulling speed and part thickness