A unified framework for parameter identification and multi-level uncertainty quantification is developed for a five-component spectral model of neuroblastoma. The model captures the interplay between aggressive and metastatic tumor phenotypes, stroma, immune response, and chemotherapy pharmacokinetics. The core novelty is the hybrid Padé-Adomian-MsDTM method, which yields analytically differentiable solutions with respect to model parameters. This enables derivation of recurrent sensitivity formulas and computation of the exact Jacobian of the truncated spectral representation. Three uncertainty quantification approaches - linearized (Delta method), Monte Carlo, and global (Sobol indices) - were implemented and compared, achieving a 5.3× speedup over classical RK45 integration. Numerical experiments on synthetic data for the aggressive Shimada morphotype confirmed reliable parameter recovery and identified tumor-infiltrating lymphocytes (TIL) and microenvironmental carrying capacity as dominant sources of prognostic uncertainty due to nonlinear threshold effects in immune dynamics. Leave-one-out cross-validation showed a predictive MAE of 7.6% on cohort-averaged clinical data. The proposed framework lays the foundation for efficient, fully differentiable personalized spectral models suitable for clinical decision support in neuroblastoma therapy.