Numerical simulations are a standard tool to investigate field theories in non-perturbative regimes. Typical algorithms used to evaluate path integrals in Euclidean space rely on importance sampling methods; i.e., a probabilistic interpretation of the Boltzmann weight eS. However, many theories of interest suffer from the infamous sign problem: the action is complex and the Boltzmann weight cannot be used as a probability distribution. Complex Langevin simulations allow numerical studies of theories that exhibit the sign problem, such as QCD at finite density. In this thesis, we study methods to investigate the phase diagram of QCD in the temperature{chemical potential plane, using the complex Langevin method. We provide results on the phase diagram for the heavy-denseapproximation of QCD (HDQCD) for three spatial volumes, using complex Langevin and the gauge cooling technique. We also present polynomial fits of the critical temperature as function of the chemical potential for each volume. Subsequently, we discuss instabilities encountered during this study, which motivated a novel technique, named Dynamic Stabilisation, which will be introduced and the theoretical ideas behind it, explained. Dynamic stabilisation was, then, used in an investigation of the dependency of the critical chemical potential on the hopping parameter. The two previous studies were used to guide a second examination of the HDQCD phase diagram, focussed around the phase boundary. Lastly, we present preliminary results on the phase diagram of QCD with fully dynamical quarks at high temperatures. This shows that complex Langevin, augmented with gauge cooling and dynamic stabilisation, is suited for investigating QCD at finite chemical potential.