In this work, soliton propagation in optical fibres is analysed by means of numerical simulation of the nonlinear Schrödinger equation, which governs optical pulse propagation in step-index monomode fibres. Since analytic solutions to this equation are admitted only for some specific cases, two numerical methods have been implemented in order to study the evolution of different kinds of input pulses, under the effects of attenuation, dispersion and nonlinearities. Showing a better performance than the Fourier Series Method in a comparative test, the Beam Propagation Method has been chosen to obtain the results here presented. Many characteristics of the fundamental, higher order and dark solitons, as well as interaction phenomena between adjacent pulses, are investigated, taking into account possible implications on optical systems performance. By properly counteracting the effects of fibre dispersion and nonlinearities, solitons can propagate without changing its shape, finding potential application in high bit-rate long distance optical communication systems.