In this article we make a local classification of n-dimensional Riemannian manifolds (M, g) with harmonic curvature and less than four Ricci eigenvalues which admit a smooth non constant solution f to the following equation $$(1)\hspace{20}{\nabla}df=f(r-{\frac{R}{n-1}}g)+x{\cdot} r+y(R)g,$$ where ∇ is the Levi-Civita connection of g, r is the Ricci tensor of g, x is a constant and y(R) a function of the scalar curvature R. Indeed, we showed that, in a neighborhood V of each point in some open dense subset of M, either (i) or (ii) below holds; (i) (V, g, f + x) is a static space and isometric to a domain in the Riemannian product of an Einstein manifold N and a static space (W, gW, f + x), where gW is a warped product metric of an interval and an Einstein manifold. (ii) (V, g) is isometric to a domain in the warped product of an interval and an Einstein manifold. For the proof we use eigenvalue analysis based on the Codazzi tensor properties of the Ricci tensor.