The reason for this mosaiked appearance is that each single image is drawn with gray scale stretched from image min to image max, and NOT across the global min/max. When you merge all images into one, well, there is only one min/max. Same is true for VRT, since VRT treats all images as one (when you have a look into the VRT file you will see the common statistics).
When both, merge and VRT, are not an option the following script might help.
In the first loop I iterate over all layers, picking the rasters and get an estimate of their min/max interval. This is the way like QGIS itself work. From these min/max values I calculate the global min/max
In the second loop the renderers for all raster layers are configured such that the gray scale is stretched across the interval global min/max.
gmin = 9999
gmax = -9999
layers = 
# loop over all layers, take rasters and estimate min, max values
for layer in iface.mapCanvas().layers():
if isinstance(layer, QgsRasterLayer):
# change percentages and sample size to increase or decrease accuracy
min_max = layer.dataProvider().cumulativeCut(1, 0.02, 0.98, theSampleSize=250000)
gmin = min(gmin, min_max)
gmax = max(gmax, min_max)
# for all rasters create a single band gray scale renderer with
# gray scale stretched across the interval [gmin..gmax]
for rasterlayer in layers:
renderer = QgsSingleBandGrayRenderer(rasterlayer.dataProvider(), 1)
# take the first band (0)
ce = QgsContrastEnhancement(rasterlayer.dataProvider().dataType(0))
# refresh canvas to show changes