This paper presents a study of fork--join systems. The fork--join system breaks down each customer into numerous tasks and processes them on separate servers. Once all tasks are finished, the customer is considered completed. This design enables the efficient handling of customers. The customers enter the system in a MAP flow. This helps create a more realistic and flexible representation of how customers arrive. It is important for modeling various real-life scenarios. Customers are divided into $K \ge 2$ tasks and assigned to different subsystems. The number of tasks matches the number of subsystems. Each subsystem has a server that processes tasks, and a buffer that temporarily stores tasks waiting to be processed. The service time of a task by the k-th server follows a PH (phase-type) distribution with an irreducible representation ($\beta_k$, $S_k$), $1 \le k \le K$. An analytical solution was derived for the case of $K=2$ when the input MAP flow and service time follow a PH distribution. We have efficient algorithms to calculate the stationary distribution and performance characteristics of the fork--join system for this case. In general cases, this paper suggests using a combination of Monte Carlo and machine learning methods to study the performance of fork--join systems. In this paper, we present the results of our numerical experiments.