According to recent experimental and numerical investigations if the characteristic size of a specimen is in the submicron size regime several new interesting phenomena emerge during the deformation of the samples. Since in such a systems the boundaries play a crucial role, to model the plastic response of submicron sized crystals it is crucial to determine the dislocation distribution near the boundaries. In this paper a phase field type of continuum theory of the time evolution of an ensemble of parallel edge dislocations with identical Burgers vectors, corresponding to the dislocation geometry near boundaries, is presented. Since the dislocation-dislocation interaction is scale free (1/r), apart from the average dislocation spacing the theory cannot contain any length scale parameter. As shown, the continuum theory suggested is able to recover the dislocation distribution near boundaries obtained by discrete dislocation dynamics simulations. One may expect, however, that for a large number of problems not all the details accounted by DDD simulations are important, the response of the dislocation network can be well described on a continuum level. Although several such continuum theories of dislocations have been developed, [27-37] most of them correspond either to mean field approximation or are based on completely phenomenological grounds. However, the role of dislocation-dislocation correlation, crucial because of the long range nature of dislocation-dislocation interaction, is far from understood. Correlation effects are taken into account in a systematic manner only in the limit when the signed dislocation density κ (geometrically necessary dislocation (GND) density) is much smaller than the stored density ρ. [38][39][40][41][42].With the advance of nanotechnology the characteristic size of the microstructure of crystalline materials reduced to the submicron level. As a consequence, the role of boundaries (sample surface, grain boundary, etc.) has become even more important than earlier. So, to model the plastic response of samples with features on the submicron scale it is crucial to determine the dislocation distribution near the boundaries. Close to a boundary the GND density is often comparable to the stored one, so the assumption |κ| ρ is not valid. The dislocation distribution near a boundary is traditionally described by the 1D pile-up of the dislocations [43]. For many real dislocation configurations, however, the interaction between dislocations in different slip planes is important requiring to go up to modeling in minimum of 2D. In this paper a phase field type theory is suggested for the simplest possible 2D dislocation arrangement consisting of straight parallel dislocations with single slip. The evolution equations of the dislocation densities are obtained from a functional of the dislocation densities and the stress potential. In contrast to other approaches suggested recently, where a set of walls of dislocations with equidistant slip distances is considered to model the dislocation config...