A search for a pair of neutral, scalar bosons with each decaying into two W bosons is presented using 36.1 fb−1 of proton-proton collision data at a centre-of-mass energy of 13 TeV recorded with the ATLAS detector at the Large Hadron Collider. This search uses three production models: non-resonant and resonant Higgs boson pair production and resonant production of a pair of heavy scalar particles. Three final states, classified by the number of leptons, are analysed: two same-sign leptons, three leptons, and four leptons. No significant excess over the expected Standard Model backgrounds is observed. An observed (expected) 95% confidence-level upper limit of 160 (120) times the Standard Model prediction of non-resonant Higgs boson pair production cross-section is set from a combined analysis of the three final states. Upper limits are set on the production cross-section times branching ratio of a heavy scalar X decaying into a Higgs boson pair in the mass range of 260 GeV ≤ mX ≤ 500 GeV and the observed (expected) limits range from 9.3 (10) pb to 2.8 (2.6) pb. Upper limits are set on the production cross-section times branching ratio of a heavy scalar X decaying into a pair of heavy scalars S for mass ranges of 280 GeV ≤ mX ≤ 340 GeV and 135 GeV ≤ mS ≤ 165 GeV and the observed (expected) limits range from 2.5 (2.5) pb to 0.16 (0.17) pb.