We present a new approach to calculate the particle distribution function about relativistic shocks including synchrotron losses using the method of lines with an explicit finite difference scheme. A steady, continuous, one dimensional plasma flow is considered to model thick (modified) shocks, leading to a calculation in three dimensions plus time, the former three being momentum, pitch angle and position. The method accurately reproduces the expected power law behaviour in momentum at the shock for upstream flow speeds ranging from 0.1c to 0.995c (1 < \Gamma < 10). It also reproduces approximate analytical results for the synchrotron cutoff shape for a non-relativistic shock, demonstrating that the loss process is accurately represented. The algorithm has been implemented as a hybrid OpenMP--MPI parallel algorithm to make efficient use of SMP cluster architectures and scales well up to many hundreds of CPUs.