A novel spectral model of neuroblastoma dynamics is proposed. This model is based on an extended nonlinear "tumor–immunity" system of equations, accounting for tumor heterogeneity, the immune response, and the pharmacokinetics of therapy. To address the stiff system, a hybrid Padé [2/2]-Adomian-MsDTM approach is employed, which reformulates the problem as recurrent algebraic relations. The model parameters are linked to morphological and immunological characteristics (Shimada classification), enabling patient-specific adaptation without increasing system dimensionality. The Padé approximation is shown numerically to expand the convergence domain and suppress boundary oscillations under sharp therapeutic interventions. The model provides a computationally efficient framework for simulating disease dynamics and evaluating individualized chemotherapy regimens; clinical validation on patient cohorts is an identified next step.