Predicting the antitumor responses of cytotoxic T-lymphocytes (i.e., CD8+T cells) can be accomplished by numerical simulations of a mathematical model of tumor-immune cell interactions to control tumor growth in cancer therapy. Cytotoxic T lymphocytes are considering as the main effector cells to target tumor cells that are in a dormant avascular state through the chemotactic mechanism and also during the cytotoxic activity of immune cells. We study a dynamical model that is formulated based on the nonlinear partial differential equations describing the interactions between tumor-infiltrating cytotoxic lymphocytes, tumor cells, TICL-tumor cell complexes, inactivated tumor-infiltrating cytotoxic lymphocytes, dead tumor cells, and chemokine. Here, a new meshless method, i.e., the direct radial basis function partition of unity (RBF-PU) technique for the spatial discretization and a semi-implicit backward differentiation formula of order 1 for the temporal discretization are employed. The proposed meshless technique is a generalized form of the original RBF-PU method such that the computational cost of computing the derivatives could be reduced in comparison with its standard form. Besides, the method depends only on the location of nodes located on the domain of the problem for the approximation, where both types of radial functions, i.e., a conditionally positive definite function, i.e., the polyharmonic splines and a positive definite function, namely Matérn function have been used The obtained fully discrete difference scheme based on the new spatial approximation is solved using an iterative method, i.e., the generalized minimal residual method with a zero-fill incomplete lower–upper preconditioner. Finally, we present some simulation results, in which we conclude that the tumor’s growth entirely depends on CTLs activation, and this behavior can be predicted by simulating the proposed mathematical model, which may have a significant role in improving cancer immunotherapy. Besides, the developed numerical method is more accurate than the classical finite element method for solving the studied model as shown in one example.