We explore the applicability of MATLAB for 3D computational fluid dynamics (CFD) of shear-driven indoor airflows. A new scale-resolving, large-eddy simulation (LES) solver titled DNSLABIB is proposed for MATLAB utilizing graphics processing units (GPUs). In DNSLABIB, the finite difference method is applied for the convection and diffusion terms while a Poisson equation solver based on the fast Fourier transform (FFT) is employed for the pressure. The immersed boundary method (IBM) for Cartesian grids is proposed to model solid walls and objects, doorways, and air ducts by binary masking of the solid/fluid domains. The solver is validated in two canonical reference cases. Then, we demonstrate the validity of DNSLABIB in a room geometry by comparing the results against another CFD software (OpenFOAM). Next, we demonstrate the solver performance in several isothermal indoor ventilation configurations and the implications of the results are discussed in the context of airborne transmission of COVID-19. The novel numerical findings using the new CFD solver are as follows. First, a linear scaling of DNSLABIB is demonstrated and a speed-up by a factor of 3-4 is also demonstrated in comparison to similar OpenFOAM simulations. Second, ventilation in three different indoor geometries are studied at both low (0.1m/s) and high (1m/s) airflow rates corresponding to Re = 5000 and Re = 50000. An analysis of the indoor CO 2 concentration is carried out as the room is emptied from stale, high CO 2 content air. We estimate the air changes per hour (ACH) values for three different room geometries and show that the numerical estimates from 3D CFD simulations may differ by 80-150 % (Re = 50000) and 75-140 % (Re = 5000) from the theoretical ACH value based on the perfect mixing assumption. Third, the analysis of the CO 2 probability distributions (PDFs) indicates a relatively non-uniform distribution of fresh air indoors. Fourth, utilizing a time-dependent Wells-Riley analysis, an example is provided on the growth of the cumulative infection risk which is shown to reduce rapidly after the ventilation is started. The average infection risk is shown to reduce by a factor of 2 for lower ventilation rates (ACH=3.4-6.3) and 10 for the the higher ventilation rates (ACH=37-64). Finally, we utilize the new solver to comment on respiratory particle transport indoors. The primary contribution of the paper is to provide an efficient, GPU compatible CFD solver environment enabling scale-resolved simulations (LES/DNS) of airflow in large indoor 1 geometries on a desktop computer. The demonstrated efficacy of MATLAB for GPU computing indicates a high potential of DNSLABIB for various future developments on airflow prediction.