A numerical method to solve the set of Dyson‐like equations in the framework of non‐equilibrium Green's functions is presented. The technique is based on the self‐consistent solution of the Dyson equations for the interacting single‐particle Green's function once a choice for the self‐energy, functional of the single‐particle Green's function itself, is made. The authors briefly review the theory of the non‐equilibrium Green's functions in order to highlight the main point useful in discussing the proposed technique. Then, the numerical implementationis presented and discussed, which is based on the distribution of the Keldysh components of the Green's function and the self‐energy on a grid of processes. It is discussed how the structure of the considered self‐energy approximations influences the distribution of the matrices in order to minimize the communication time among processes and which should be considered in the case of other approximations. The authors give an example of the application of this technique to the case of quenches in ultracold gases and to the single impurity Anderson model, also discussing the convergence and the stability features of the approach.