This study aims at extending and implementing an iterative spectral scheme for fractional-order unsteady nonlinear integro-partial differential equations with weakly singular kernel. In this scheme, the unknown function u(x, t) is estimated by using shifted Gegenbauer polynomials vector Λ(x, t), and Picard iterative scheme is used to handle underlying nonlinearities. Some novel operational matrices are developed for the first time in order to approximate the singular integral like, ${\int }_0^x{\int }_0^{y}u(p{a}_1+b_1,q{a}_2+b_2,t)/{(x^{{\rho }_1}-p^{{\rho }_1})}^{{\alpha }_1}{(y^{{\rho }_2}-q^{{\rho }_2})}^{{\alpha }_2}\text{d}q\text{d}p$ \end{document} and ${\int }_0^t{u}^{\gamma }(\mathbf{x},\xi )/{(t^{{\rho }_3}-{\xi }^{{\rho }_3})}^{{\alpha }_3}\text{d}\xi$, where ρ's > 1, 0 < α's < 1 by means of shifted Gegenbauer polynomials vector. The advantage of this extended method is its ability to convert nonlinear problems into systems of linear algebraic equations. A computer program in Maple for the proposed scheme is developed for a sample problem, and we validate it to compare the results with existing results. Six new problems are also solved to illustrate the effectiveness of this extended computational method. A number of simulations are performed for different ranges of the nonlinearity n, α, fractional-order, ρ, and convergence control M, parameters. Our results demonstrate that the extended scheme is stable, accurate, and appropriate to find solutions of complex problems with inherent nonlinearities.