## 4 Results and Discussion  
The bvp4c boundary value problem solver using MATLAB software was used to solve the governing nonlinear ordinary differential equation for boundary layer flow and heat transfer with its boundary condition. The three-stage Lobatto IIIa method is a numerical solution in bvp4c solver. The Lobatto IIIa method was proposed by [PERSON]. This method provides specified accuracy, where an initial guess is supplied at initial mesh points and changes step-size. According to [PERSON] [21], the bvp4c solver produced satisfactory solutions compared with the shooting method and Keller box. The obtained results were validated by comparing the output with previously reported results, and it was shown that the comparison was in an excellent agreement with the published results. These numerical results are presented in the form of tables and graphs. Next, parameter values \(\varphi_{2}\), \(M\), \(\epsilon\), \(Rd\), and nanofluid types were selected to obtain numerical results for variation \(\lambda\) in shrinking case \(\lambda<0\). In this study, the Prandtl number value was set as \(\text{Pr}=6.2\), which shows that the water and nanoparticle types selected were carbon nanotubes (CNTs), which were single-wall CNTs and multi-wall CNTs, as well as copper (Cu). Next, the \(S\) parameter shows that the model is permeable and is considered to improve the suction rate, \(S>0\), and that the second solution is obtained when in a certain range, which is \(S>2\). Based on this study, this model depends on suction strength to control boundary layer separation on molecules that move slowly in flow. Meanwhile, the various numerical results of \(f^{\prime\prime}(0)\) and \(-\theta^{\prime}(0)\) were compared with previous results reported by the authors of the related studies and are shown in Table 3. Based on Table 3, each numerical result of the first and second solutions were compared and are shown to have achieved a good agreement, meaning that the numerical values of the present results are almost similar to the values of previous results, which are from [29; 30; 31]. The precision of the solutions to the boundary value problems found by using MATLAB bvp4c is comparable with previous studies, proving the compatibility of the numerical results in Table 3.  
Apart from that, Table 4 shows numerical results for \(-\theta^{\prime}(0)\), which are compared with previous results [31; 32; 33; 34; 35]. Based on this table, it is proven that the present numerical results are almost identical to previous results and that they achieve good agreement. Based on the present results in Table 4, the higher value of the Prandtl number shows the increasing \(-\theta^{\prime}(0)\) values, which causes the temperature boundary layer thickness to become thinner, demonstrating the improvement in the heat transfer process. Thus, heat transfer enhancement becomes better. In general, Figures 2-5 show the existence of the dual solution clearly, where the randomly shrinking case, \(\lambda_{c}\leq\lambda\leq-\lambda\), can be observed. The unique solution exists at point \(\lambda_{c}=\lambda\), where \(\lambda_{c}\) is the critical point. In region \(\lambda_{c}\leq\lambda\leq-\lambda\), the Tiwari-Das model was developed to describe the behavior of the MHD hybrid nanofluid flow and heat transfer on the permeable exponentially shrinking surface with a few parameters. A numerical solution does not exist in region \(\lambda_{c}>\lambda\) because this case describes the incompatibility of the mathematical model in this region. The boundary condition separation from the exponentially uneven surface shows that the shrinking case causes the Tiwari-Das model, which was developed and used based on boundary layer approximation theory, to be invalid. The second solution obtained is limited over a certain range of \(\lambda\) values, which is because the specific value of \(\lambda\) cannot track the exact solution according to the relative error