<p>The numerical solution of time-dependent systems is often carried out using a step-wise updating procedure known as marching in time or simply time-marching methods. They require completing previous steps to yield an update for the next time step. This imposes a barrier in implementing the algorithm in a computer program, which prevents parallel speed-up by increasing the serial fraction of the program. This paper analyzes a method based on the space-time discretization of a time-dependent system using Discrete Picard Iterations (DPI). Unlike time marching schemes, the updating barrier is removed in the current formulation by converting the solution procedure into a set of larger matrix–vector products that correspond to the simultaneous computation of spatial residuals at various time levels. The stability and accuracy of the proposed method are investigated. The accuracy and computation time are compared to various popular time-marching methods involving structured and higher order/spectral unstructured grids in one and two-dimensional space. It is shown that when the inverse of the Chebyshev differentiation matrix is used in the DPI framework, the resulting scheme is superior in both accuracy and CPU time. Various forms of iterative methods and preconditionings are implemented and evaluated. This includes GMRES, IDR(s), and Conjugate Gradient methods. It is shown that the stabilized Bi-conjugate gradient method produces the best iterative solver when a band-LU preconditioner is used. This method, which relies heavily on matrix multiplication, is strategically consistent with the recent development in the hardware implementation of matrix multiplication accelerators on ASICs and FPGAs.</p>