In this thesis, a numerical study of double layers in plasma is presented. Double layersare structures consisting of two oppositely charged space charge layers, creating a finitechange in electrostatic potential over the double layer. The motivation for the topic isfound in the general interest for magnetospheric physics, where double layers are knownto accelerate particles into the ionosphere, creating the aurora, as seen at higher latitudes. In the first part of this thesis, theories and observations of double layers as a plasmaphysics phenomena are studied. Two double layer related phenomena, the Bunemaninstability and adiabatic cooling, are presented and their expected influence on doublelayers given. Then, through the introduction of BGK-solutions and the Water-bag model,a one dimensional, time stationary, electrostatic model of double layers is formed. In the second part of this thesis, the theories discussed in the first part are imple-mented into a one dimensional numerical model of double layers. The numerical modelemphasizes strong (i.e. Buneman regime) double layers, and their existence is investi-gated for a large number of different plasma and boundary conditions. The results of the numerical model can confirm important observations made by previ-ous numerical and experimental studies, like the scaling law, Bohm criterion and presenceof two-stream (i.e. Buneman) instabilities. The present numerical model improves onthese studies by introducing simulations with varying combinations of boundary con-ditions, and a wider selection of plasma conditions for which double layers have beensimulated. For the numerical studies of this thesis, a one dimensional Particle-in-Cell code wasdesigned and developed. A thorough description of the program is given in the main partof this thesis, and the source code is shown in the appendix.