We present BILL2D, a modern and efficient C++ package for classical simulations of two-dimensional Hamiltonian systems. BILL2D can be used for various billiard and diffusion problems with one or more charged particles with interactions, different external potentials, an external magnetic field, periodic and open boundaries, etc. The software package can also calculate many key quantities in complex systems such as Poincare sections, survival probabilities, and diffusion coefficients. While aiming at a large class of applicable systems, the code also strives for ease-of-use, efficiency, and modularity for the implementation of additional features. The package comes along with a user guide, a developer's manual, and a documentation of the application program interface (API). Program summary Program title: Bill2d Catalogue identifier: AEYLv1_0 Program summary URL: http://cpc.cs.qub.ac.uk/summaries/AEYLv1_0.html Program obtainable from: CPC Program Library, Queen's University, Belfast, N. Ireland Licensing provisions: GNU General Public License, version 3 No. of lines in distributed program, including test data, etc.: 37098 No. of bytes in distributed program, including test data, etc.: 1155037 Distribution format: tar.gz Programming language: C++(14). Computer: Tested on x86 and x86 64 architectures. Operating systems: Tested on Linux, and OS X versions 10.9-10.11. Has the code been vectorized or parallelized?: Shared memory parallelization when simulating ensembles of systems. Vectorization of operations with R-2 vectors. RAM: Simulation dependent: kilobytes to gigabytes Classification: 4.3, 7.8, 7.9, 7.10, 16.9. External routines: Boost, CMake, GSL, HDF5; and optionally Google-Mock, GoogleTest, and Doxygen Nature of problem: Numerical propagation of classical two-dimensional single and many-body systems, possibly in a magnetic field, and calculation of relevant quantities such as Poincare sections, survival probabilities, diffusion co-efficients, etc. Solution method: Symplectic numerical integration of Hamilton's equations of motion in Cartesian coordinates, or solution of Newton's equations of motion if in a magnetic field. The program implements several well-established algorithms. Restrictions: Pointlike particles with equal masses and charges, although the latter restrictions are easy to lift. Unusual features: Program is efficient, extremely modular and easy to extend, and allows arbitrary particle-particle interactions. Additional comments: The source code is also available at https://bitbucicet.orgisolanpaa/bill2d. See README for locations of user guide, developer manual, and API docs. Running time: From milliseconds to days, depends on type of simulation. (C) 2015 Elsevier B.V. All rights reserved.