A fully spectral multi-domain method has been developed and applied to three applications within ideal MHD, compressible Navier-Stokes, and a two-fluid plasma turbulence model named the Weiland model. The time-spectral method employed is the Generalized Weighted Residual Method (GWRM), where all domains such as space, time, and parameter space are spectrally decomposed with Chebyshev polynomials. The spectral decomposition of the temporal domain allows the GWRM to reach spectral accuracy in all dimensions. The GWRM linear/nonlinear algebraic equations are solved using an Anderson Acceleration (AA) method and a newly developed Quasi Semi-Implicit root solver (Q-SIR). Up to 85% improved convergence rate was obtained for Q-SIR as compared to AA and in certain cases only Q-SIR converged. In the most challenging simulations, featuring steep gradients, the GWRM converged for time intervals roughly two times larger than typical time steps for explicit time-marching schemes, being limited by the CFL condition. Time intervals up to 70 times larger than those of explicit time-marching schemes were used in smooth regions. Furthermore, the most computationally expensive algorithm, namely the product of two Chebyshev series, has been GPU accelerated with speedup gains of several thousands compared to a CPU.