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\ \mathrm{\mathcal{B}}\left({B}_s^0\to {\mu}^{+}{\mu}^{-}\right)=\left({3.2}_{-1.0}^{+1.1}\right)\times {10}^{-9} \ and an upper limit ℬB0→μ+μ−<4.3×10−10\ \mathrm{\mathcal{B}}\left({B}^0\to {\mu}^{+}{\mu}^{-}\right)<4.3\times {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\ \mathrm{\mathcal{B}}\left({B}_s^0\to {\mu}^{+}{\mu}^{-}\right)=\left({2.8}_{-0.7}^{+0.8}\right)\times {10}^{-9} \ and ℬB0→μ+μ−<2.1×10−10\ \mathrm{\mathcal{B}}\left({B}^0\to {\mu}^{+}{\mu}^{-}\right)<2.1\times {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→μ+μ−\ \mathrm{\mathcal{B}}\left({B}^0\to {\mu}^{+}{\mu}^{-}\right)-\mathrm{\mathcal{B}}\left({B}_s^0\to {\mu}^{+}{\mu}^{-}\right) \ plane.