The motion of a passive spherical particle in a fluid has been widely described via a balance of force equations known as a Generalized Langevin Equation (GLE) where the covariance of the thermal force is related to the time memory function of the fluid. For viscous fluids, this relationship is simply a delta function in time, while for a viscoelastic fluid it depends on the constitutive equation of the fluid memory function. In this paper, we consider a general setting for linear viscoelasticity which includes both solvent and polymeric contributions, and a family of memory functions known as the generalized Rouse kernel. We present a statistically exact algorithm to generate paths which allows for arbitrary large time steps and which relies on the numerical evaluation of the covariance of the velocity process. As a consequence of the viscoelastic properties of the fluid, the particle exhibits subdiffusive behavior, which we verify as a function of the free parameters in the generalized Rouse kernel. We then numerically compute the mean first passage time of a passive particle through layers of different widths and establish that, for the generalized Rouse kernel, the mean first passage time grows quadratically with the layer’s width independently of the free parameters. Along the way, we also find the linear scaling of the mean first passage time for a layer of fixed width as a function of the particle’s radius.