This paper presents a mathematical analysis of bioconvection heat, mass, and motile microorganisms transfer over a stretching sheet in a medium filled with a fluid containing gyrotactic microorganisms. Cross-diffusion is taken into account in the medium. Using the boundary layer approximations, the set of unsteady partial differential equations governing the fluid flow is transformed into nonlinear PDEs form and then solved numerically using bivariate spectral relaxation method (BSRM) and bivariate spectral quasi-linearization method (BSQLM). A comparison between BSRM and BSQLM is made for the first time in this work. The accuracy and convergence analysis of the methods are also discussed. The methods are found to be convergent and give very accurate results with very few grid points in the numerical discretization procedure. A parametric study of the entire flow regime is carried out to illustrate the effects of various governing parameters on the fluid properties and flow characteristics. The results obtained show a significant effect of cross-diffusion on the fluid properties and flow characteristics. The Dufour number was found to increase the local Sherwood number and density number of motile microorganisms while decreasing the Nusselt number, and the reverse effect is true for the Soret number. Furthermore, the Nusselt number, Sherwood number, and density number of motile microorganisms are highly influenced by buoyancy and bioconvection parameters.