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,\tag1$$ 
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 curvature, $g$ has harmonic curvature.  (This is a well-known special case of harmonic curvature.  In dimension $3$, harmonic curvature is *equivalent* to being conformally flat with constant scalar 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$. For example, when $n=2$, the equation takes the form
$$
u(t)\,f''(t)-2\,u(t)u'(t)\,f'(t)+2\,u'(t)^2-u(t)u''(t)=0.
\tag2
$$  
Thus, $f$ can be found by quadrature once the function $u$ has been specified.  For the generic pair $(u,f)$ satisfying (1) and (2), the corresponding $\lambda(t) = \bigl(u'(t)f'(t)-u''(t)\bigr)u(t)^{-3}$ will be 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 fact, when $k>0$, there are periodic solutions $u(t)$ defined on the entire line $I=\mathbb{R}$ to the equation corresponding to (1), and one can use this to construct examplees defined on compact manifolds.)

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.