We present a systematic approach to the simulation of EPR spectra of several different spin systems with spin S > 1/2 at arbitrary microwave frequencies and various zero-field splittings (zfs) using direct diagonalization of the complete spin energy matrix. With this method effective simulations can be obtained for systems at every microwave energy irrespective of the zfs interactions and it is therefore also expected to work for systems in the ‘low field limit’ at low microwave frequencies when the zfs energy dominates and where perturbation approaches fail. Here we have used metal complexes in biological systems with substantial zfs interactions, especially the nitrogenase ‘super-clusters’, in order to test the simulation routine. The nitrogenase MoFe-protein contains two types of very complex clusters, FeMo-cofactor and the P-cluster, that can be prepared in several different paramagnetic, oxidation states. It is therefore an ideal system for probing the usefulness of the method. Good spectral simulations were obtained for the FeMo-cofactor signal, with total spin S = 3/2, and the P-cluster signals, with total spins S = 7/2 and S = 9/2 and probably S = 5/2, from the nitrogenase MoFe-protein of Klebsiella pneumoniae (Kp1) at various oxidation levels. The simulations of the S = 3/2 FeMo-cofactor signal reproduces both the ΔM_S = 1 transitions generating a strongly rhombic signal from transitions between the M_S = ±1/2 spin levels and the ΔM_S = 3 transition within the M_S = ±3/2 Kramers doublet that gives rise to a signal at g = 6.0. Ferricyanide-oxidised Kp1 gives us the rare opportunity to observe signals from both S = 7/2 and an earlier putatively assigned S = 9/2 spin state from a cluster in a biological environment. In contrast to earlier measurements on the thionine oxidised protein, the experimentally measured S = 7/2 and S = 9/2 states are observed simultaneously for ferricyanide oxidised Kp1 and the individual spectra have been simulated.