A method of calculating the time correlation functions for electron spin is proposed, based on solving the time-dependent Schrodinger equation for a spin Hamiltonian that contains a term varying randomly in discrete time steps. It is applied to the study of electron spin relaxation in aqueous solution of nickel(II) ions with S=1. The random term in the spin Hamiltonian in this case is the zero-field splining (ZFS) interaction. The method is evaluated by an application to a model system (the pseudorotation model) for which an analytical solution to the electron spin relaxation problem is known. The same method is then employed to study the electron and nuclear spin dynamics in a system where the time variation of the zero-field splitting is obtained by a combination of ab initio quantum chemistry and molecular dynamics simulations.