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, 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. © 2012 The Authors Monthly Notices of the Royal Astronomical Society © 2012 RAS.