A semi-Lagrangian method is developed for the solution of Lagrangian transport equations and stochastic differential equations that is consistent with Discontinuous Spectral Element Method (DSEM) approximations of Eulerian conservation laws. The method extends the favorable properties of DSEM that include its high-order accuracy, its local and boundary fitted properties and its high performance on parallel platforms for the concurrent semi-Lagrangian and Eulerian solution of a class of time-dependent problems that can be described by coupled Eulerian-Lagrangian formulations. Such formulations include the probabilistic models used for the simulation of chemically reacting turbulent flows or particle-laden flows. Motivated by the high-fidelity simulation of chemically, reacting, turbulent shear layers in propulsion systems, this thesis also reports on the stability of the interaction of two cold flow shear layers. Consistent with an explicit, DSEM discretization, the semi-Lagrangian method seeds particles at Gauss quadrature collocation nodes within a spectral element. The particles are integrated explicitly in time according to a governing deterministic and/or stochastic differential equation and form the nodal basis for an advected interpolant. This interpolant is mapped back in a semi-Lagrangian fashion to the Gauss quadrature points through a least squares fit using constraints for element boundary values and optional constraints for mass and energy preservation. Stochastic samples are averaged on the quadrature nodes. The stable explicit time step of the DSEM solver is sufficiently small to prevent particles from leaving the element’s bounds. The semi-Lagrangian method is hence local and parallel and does not have the grid complexity, and parallelization challenges of the commonly used Lagrangian particle solvers in particle-mesh methods for solution of Eulerian-Lagrangian formulations. Numerical tests in one and two dimensions for advection and (stochastic) diffusion problems show that the method converges exponentially for constant and non-constant advection and diffusion velocities. The use of mass and energy constraints can improve accuracy depending on the order of accuracy of the time integrator. The linear and non-linear growth of instabilities in a two shear layer cold flow configuration is studied using a combination of Linear Stability Analysis (LSA) and Direct Numerical Simulation (DNS) based on a DSEM approximation of the first-principle Navier-Stokes equations. The linear growth of spatial and temporal LSA modes for a single shear layer is verified using DNS. Then DNS is used to study the transition from linear growth to non-linear instabilities. In DNS of a non-parallel flow with a spatially growing viscous LSA mode an absolute instability is observed. In that case, a match between DNS and LSA is difficult because of spurious modes introduced by the reflections at the boundaries. In a temporal analysis, it is found that the growth of the temporal mode is greater with an increasing difference in the strength of the two shear layers for both inviscid and viscous flows. Longer DNS runs show that presence of a stronger shear layer enhances vortex shedding and vortex pairing mechanism of a shear layer. This enhanced mixing in the non-linear region shows a correlation with the growth of perturbation in the linear region.