A study of the decays Bs0 → μ+μ− and B0 → μ+μ− has been performed using 26.3 fb−1 of 13 TeV LHC proton-proton collision data collected with the ATLAS detector in 2015 and 2016. Since the detector resolution in μ+μ− invariant mass is comparable to the Bs0-B0 mass difference, a single fit determines the signal yields for both decay modes. This results in a measurement of the branching fraction ℬBs0→μ+μ−=3.2−1.0+1.1×10−9 and an upper limit ℬB0→μ+μ−<4.3×10−10 at 95% confidence level. The result is combined with the Run 1 ATLAS result, yielding ℬBs0→μ+μ−=2.8−0.7+0.8×10−9 and ℬB0→μ+μ−<2.1×10−10 at 95% confidence level. The combined result is consistent with the Standard Model prediction within 2.4 standard deviations in the ℬB0→μ+μ−−ℬBs0→μ+μ− plane.