My answer will be related to solid (crystalline) substances. In this case you need a well oriented sample and a set of its spectra, taken at different orientations with respect to the external magnetic field.  The data you need to extract from these spectra are essentially the values of resonance fields as a function(s) of sample orientations.  When plotted in polar coordinates, they will reveal the symmetry of your copper complex environment.  Now you are ready to construct the appropriate spin hamiltonian relevant to your experiment. Once you have it, you can simulate the positions of ESR lines, again as function(s) of sample orientation.  Yes, you have to guess the values of appropriate constants, say components of g-tensor and A-tensor.

After your simulations are overlapping with experimental data - you are done.  Of course, thing could be much simpler, if the simulations could be performed in analytical way, not numerically, what is only possible in rare cases (Cu+ ion, with effective spin 1/2 may be one of them).  The rare cases are, in principle, those when you can solve the eigen equation analytically, that is when S

Cu2+ ions are normally fairly simple to model provided they are isolated and act as S = 1/2 spins. [If they come close together or have bridging ligands then magnetic exchange can make their spectra more unusual; copper acetate, Cu2(OAc)4.2H2O is a good example]. If your sample was measured in solution at room temperature you most likely have a 1:1:1:1 quartet due to coupling to the Cu2+ ion (both isotopes have I = 3/2). You may have additional hyperfine coupling to ligand atoms such as N but this is typically not well resolved. The separation of the peaks (in Gauss) corresponds to the hyperfine coupling a(Cu) and you can determine the isotropic g-value from g = (h x nu)/(beta x H) where nu is the microwave frequency, H is the applied field and beta the Bohr magneton. Freeware programs such as WINSIM will give you a reasonable estimate of both the coupling constant and g-value.

If you have a solid state sample or frozen solution then life is more complex as there is no molecular tumbling and you observe an EPR spectrum which is the sum of all orientations rather than the average. In this case we resolve the g-tensor into three principle components g(x), g(y) and g(z). For Cu2+ we often find that a Jahn-Teller distortion provides a pseudo-axial spectrum in which g(x) ~g(y) g(z). In addition the hyperfine coupling is also split into three components A(Cu,x), A(Cu,y) and A(Cu,z). Technically the g and A-matrices need not be coincident but, for many systems they are sufficiently close that you can get a reasonable simulation assuming they are coincident. The issue is probably to find a simulation program which will simulate anisotropic spectra. If you have access to Simfonia (Bruker software I think) then this will cope but I don't personally like the interface.There is some other software which I recall is called easyspin but uses matlab (which I don't have). I personally use an old fortran program called PIP written by M. Nilges at the Illinois EPR centre for which I wrote a windows interface some years ago. If you have access to windows XP(!) or a windows XP emulator in windows 7, I can send you the code and instructions to install the s/w if that would be helpful.