<p>This paper presents a novel spatiotemporal partial differential equation (PDE) model that captures the complex interactions among tumor cells, immune responses, and chemotherapy effects. The proposed framework incorporates nonlinear diffusion, saturation kinetics, and spatially distributed drug administration to realistically describe tumor immune therapy dynamics. For the numerical treatment, the system is discretized using the Legendre spectral collocation method, which provides high-order accuracy and rapid convergence with relatively low computational cost. Numerical simulations reveal that spatial heterogeneity plays a crucial role in tumor progression and treatment efficacy, significantly influencing immune infiltration and drug penetration. The results further demonstrate conditions under which complete tumor eradication, tumor relapse, or persistent tumor states may occur, depending on immune efficiency and chemotherapy intensity. Overall, the study highlights the effectiveness of spectral methods for biomedical modeling and offers a flexible mathematical foundation for the design and optimization of personalized cancer treatment strategies.</p>