In this work, a methodology for parallel numerical solution of general systems of non-linear differential and algebraic equations (ODE or DAE) on distributed memory systems is proposed and implemented. The methodology consists of the following parts: (1) an algorithm for transformation of model equations into bytecode instructions suitable for parallel evaluation on diverse types of computing devices, (2) data structures for model specification, (3) an algorithm for partitioning of general systems of equations, (4) an algorithm for inter-process data exchange, and (5) simulation software for integration of general DAE systems in time. The model specification contains only the low-level model description with the minimum information required for integration in time. This way, the model specification data structures are used as a simple platform-independent binary interface for model exchange. The partitioning algorithm accepts four different balancing constraints that can be used to precisely balance the computation and memory loads in critical phases of the numerical solution. For computationally intensive tasks the software utilises multi-level parallelism techniques such as hybrid MPI/OpenMP and heterogeneous MPI/OpenCL. The proposed methodology is applied to a medium-scale transient two-dimensional phase separation process. Six different phases of the numerical solution are analysed. Nine different strategies for load balancing are applied. Simulation results, an overall performance and performance of individual phases, an efficiency of the preconditioner, the quality of the load balancing prediction, and overheads due to the load imbalance are discussed in details.