We consider the type I multiple orthogonal polynomials (MOPs) (An,m, Bn,m), deg An,m ≤ n − 1, deg Bn,m ≤ m − 1, and type II MOPs Pn,m, deg Pn,m = n + m, satisfying non-hermitian orthogonality with respect to the weight e −z 3 on two unbounded contours γ 1 and γ 2 on C, with (in the case of type II MOPs) n conditions on γ 1 and m on γ 2 . Under the assumption that n, m → ∞, n n + m → α ∈ (0, 1)we find the detailed (rescaled) asymptotics of An,m, Bn,m and Pn,m on C, and describe the phase transitions of this limit behavior as a function of α. This description is given in terms of the vector critical measure µα = (µ 1 , µ 2 , µ 3 ), the saddle point of an energy functional comprising both attracting and repelling forces. This critical measure is characterized by a cubic equation (spectral curve), and its components µ j live on trajectories of a canonical quadratic differential on the Riemann surface of this equation. The structure of these trajectories and their deformations as functions of α was object of study in our previous paper [Adv. Math. 302 (2016), 1137-1232, and some of the present results strongly rely on the analysis carried out there.We conclude that the asymptotic zero distribution of the polynomials An,m and Pn,m are given by appropriate combinations of the components µ j of the vector critical measure µα. However, in the case of the zeros of Bn,m the behavior is totally different, and can be described in terms of the balayage of the signed measure µ 2 − µ 3 onto certain curves on the plane. These curves are constructed with the aid of the canonical quadratic differential , and their topology has three very distinct characters, depending on the value of α, and are obtained from the critical graph of .Once the trajectories and vector critical measures are studied, the main asymptotic technical tool is the Deift-Zhou nonlinear steepest descent analysis of a 3 × 3 matrix-valued Riemann-Hilbert problem characterizing both (An,m, Bn,m) and Pn,m, which allows us to obtain the limit behavior of both types of MOPs simultaneously.We illustrate our findings with results of several numerical experiments, explain the computational methodology, and formulate some conjectures and empirical observations based on these experiments.