In this thesis we develop techniques to efficiently solve numerical Partial Differential Equations (PDEs) using Graphical Processing Units (GPUs). Focus is put on both performance and re--usability of the methods developed, to this end a library, cuSten, for applying finite--difference stencils to numerical grids is presented herein. On top of this various batched tridiagonal and pentadiagonal matrix solvers are discussed. These have been benchmarked against the current state of the art and shown to improve performance in the solution of numerical PDEs. A variety of other benchmarks and use cases for the GPU methodologies are presented using the Cahn--Hilliard equation as a core example, but it is emphasised the methods are completely general. Finally through the application of the GPU methodologies to the Cahn--Hilliard equation new results are presented on the growth rates of the coarsened domains. In particular a statistical model is built up using batches of simulations run on GPUs from which the growth rates are extracted, it is shown that in a finite domain that the traditionally presented results of 1/3 scaling is in fact a distribution around this value. This result is discussed in conjunction with modelling via a stochastic PDE and sheds new light on the behaviour of the Cahn--Hilliard equation in finite domains.