We analyze a Markovian queueing system with a two-stage service process, where both stages are mandatory for each customer to complete the service. The matrix method and Laplace-Stieltjes transform (LST) are used to determine the transient probabilities of the system states and to compute the performance indices of the model. A numerical study is performed to evaluate the system’s sensitivity. To validate the findings, we compare those obtained using the adaptive neuro-fuzzy inference system (ANFIS) with those derived from the matrix method. Furthermore, we formulate a nonlinear cost function with decision variables representing the service rates of both stage 1 and stage 2. This cost function is then optimized using the particle swarm optimization (PSO) algorithm to find the optimal service rates. The optimization process helps identify the most cost-efficient configuration, improving system performance and efficiency. The proposed model is applicable in various domains, such as call centers, healthcare facilities, and banks.