A search for a pair of neutral, scalar bosons with each decaying into two $W$ bosons is presented using 36.1 fb$\mathrm{<mrow\; data-mjx-texclass="ORD">-1</mrow>}$ 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 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 $\le mX\le$ 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 $\le mX\le$ 340 GeV and 135 GeV $\le mS\le$ 165 GeV and the observed (expected) limits range from 2.5 (2.5) pb to 0.16 (0.17) pb.