Here's a simple set of examples that are special cases of the family of solutions that I was referring to.  (Describing all of the examples would take longer than I have at the moment.)

Let $I\subset\mathbb{R}$ be an open interval and consider metrics on $I\times\mathbb{R}^n$ of the form 
$$
g = u(t)^2\bigl((\mathrm{d}t)^2 + (\mathrm{d}x^1)^2 + \cdots +(\mathrm{d}x^n)^2\bigr),
$$
where $u:I\to\mathbb{R}^+$ is a positive function on $I$.   The condition that $g$ have constant scalar curvature is an autonomous nonlinear second order ODE for $u$ that is solvable in terms of Jacobian elliptic functions (i.e.,  inversion of an elliptic integral).  (I won't go into details about this until I have more time.  When $n=2$, for example, the equation is $2\,u(t)\,u''(t)-u'(t)^2+3c\,u(t)^4=0$, where $6c$ is the scalar curvature of $g$.)  

For the generic solution $u$, the metric $g$ is not Einstein.  In fact, $\mathrm{Ric}(g) = a(t)\,g + b(t)\,\mathrm{d}t^2$ for functions $a$ and $b$ determined in terms of $u$.  Generically, $b$ will not be zero.  However, since $g$ is conformally flat and has constant scalar curvatue, $g$ has harmonic curvature.

Now consider a function $f:I\to\mathbb{R}$ regarded as a function on $I\times\mathbb{R}^n$ and ask that $\mathrm{Ric}(g)+\mathrm{Hess}_g(f) = \lambda(t)\,g$ for some $\lambda(t)$.  Since $\mathrm{Hess}_g(f) = p(t)\,g + q(t)\,\mathrm{d}t^2$ for some functions $p$ and $q$ determined in terms of $f$ and $u$, this is a single linear second order equation for $f$ involving $u(t)$ that always has solutions on the interval $I$ of definition of $u$.  In fact, $f$ can easily found by quadrature once the function $u$ has been specified.  The generic solution pair $(u,f)$ will have $\lambda(t)$ nonconstant, thus producing examples of the desired type.

This can be generalized slightly by replacing $\mathbb{R}^n$ and its flat metric by $(N^n,h)$ where $h$ is a metric on $N$ of constant sectional curvature $k$. The equations for $u(t)$ and $f(t)$ change slightly, but they are still solvable using elliptic functions.

In dimension $3$, it turns out that *every* conformally flat metric $(M^3,g)$ with constant scalar curvature that admits a 'Ricci potential', i.e., a function $f$ such that $\mathrm{Ric}(g)+\mathrm{Hess}_g(f) = \lambda\,g$ for some function $\lambda$ is locally of the above form for some $(N^2,h)$.  It was after I worked that out using an exterior differential analysis that I realized that the above construction would work in any dimension.