In this study, a computational matrix-collocation method based on the Lagrange interpolation polynomial is specifically streamlined to treat the fractional nonlinear terminal value problems with multiple delays, such as the Hutchinson, the Wazewska-Czyzewska and the Lasota models in biomathematics. To do this, the robust nonlinear terms of which are smoothed to be deployed in the method. The uniqueness analysis of the solution is discussed in terms of the Banach contraction principle. An error analysis technique is non-linearly theorized and applied to improve the solutions. A programme for the method is especially developed. Thus, the outcomes of five fractional model problems constrained by terminal conditions are numerically and graphically evaluated in tables and figures. Based on the investigation of the results, one can claim that the method presents a sustainable and effective mathematical procedure for the aforementioned problems.