The paper studies the classical fork-join queueing system with M|M|1 subsystems. The analysis of this system is still relevant due to the lack of exact solutions for assessing its performance characteristics if the number of subsystems exceeds two. In addition, the fork-join system is a mathematical model of parallel or distributed computing systems that have become widespread as one of the most effective methods for processing Big Data. An approach based on graphical analysis, non-linear regression, and the use of the Nelder-Mead optimization method is proposed to estimate the mathematical expectation and dispersion of the response time of a fork-join system. As a result, the authors managed to modify the known approximations and significantly (many times) improve their approximation quality. The paper also examines the quality of the experimental data of simulation modeling used to estimate the approximation error of the obtained expressions. As a rule, this issue remains outside the scope of ongoing research in the field of this topic due to the complexity of such an analysis. And sometimes, it is due to the underestimation of the importance of this issue. The article proposes an approach to finding confidence intervals for simulation results. It provides an algorithm for their construction and also gives some recommendations.