We consider the numerical solution of a phase field model for polycrystallization in the solidification of binary mixtures in a domain Ω ⊂ R 2 . The model is based on a free energy in terms of three order parameters: the local orientation Θ of the crystals, the local crystallinity φ, and the concentration c of one of the components of the binary mixture. The equations of motion are given by an initial-boundary value problem for a coupled system of partial differential equations consisting of a regularized second order total variation flow in Θ, an L 2 gradient flow in φ, and a W 1,2 (Ω) * gradient flow in c. Based on an implicit discretization in time by the backward Euler scheme, we suggest a splitting method such that the three semidiscretized equations can be solved separately and prove existence of a solution. As far as the discretization in space is concerned, the fourth order Cahn-Hilliard type equation in c is taken care of by a C 0 Interior Penalty Discontinuous Galerkin approximation which has the advantage that the same finite element space can be used as well for the spatial discretization of the equations in Θ and φ. The fully discretized equations represent parameter dependent nonlinear algebraic systems with the discrete time as a parameter. They are solved by a predictor corrector continuation strategy featuring an adaptive choice of the time-step. Numerical results illustrate the performance of the suggested numerical method.