I am trying to solve the steady-state solution of the 1D heat equation with a known source/sink term and non-homogeneous Neumann boundary conditions, however I am not sure if an analytical solution does exist.
To summarize the equation is: $\kappa T_{zz} + f(z)=0$
and the boundary condition (BC) at the surface is: $\kappa T_z(z=0)=q_0$
For coherence, as the objective is to find the steady state solution, I assumed that the heat flux at the bottom ($q_h$) should balance the heat flux at the surface ($q_0$) plus the source term ($f(z)$) (i.e., equilibrium conditions), so that:
$\int_0^h \kappa T_{zz}\, dz= q_h - q_0 = - \int_0^h f(z)dz\,,$
the bottom boundary condition is then taken as: $\kappa T_z(z=h)=q_h=q_0- \int_0^h f(z)dz$
So, the first question is if this approach is correct.
Then, with these Neumann boundary conditions, integrating the equation one time one obtains:
$\kappa T_z = - \int f(z)dz + c_1$
from which $c_1$ can be calculated e.g., applying the boundary condition at the surface as $c_1= q_0 + \int f(z)dz\bigg|_{z=0}$.
Now, I do not know how to proceed as integrating the equation one more time will introduce the second constant $c_2$ that I do not know how to determine. So the second question is: is there an analytical solution for this problem? According to this previous question I would say that a solution exists but "the solution is not unique without imposing additional constraints". Here the constraints that (I think) one should use is the definition of $q_h$, but I am not sure if this is correct and in case how to use it.
About your first question, existence of a solution to the Neumann problem requires the total flux to be zero. To help with keeping signs straight, note that the underlying heat equation you're looking at is $\frac{\partial u}{\partial t}=\kappa \Delta u + f$, so the source is $f$. So the inward flux will be $\kappa T_z(h)-\kappa T_z(0)+\int_0^h f(z) dz$. That must be zero in order for a solution to exist (note the inclusion of $\kappa$'s). Physically, when this condition does not hold, $\int_0^h T(t,z) dz$ changes (where $T(t,z)$ is the dynamical solution to the corresponding heat equation). This precludes the existence of any steady state.
In terms of finding the solution, it is just an ordinary differential equation BVP. You can integrate it twice and plug in one boundary condition to get one of the constants of integration. The other one cannot be determined, because the solution to the problem is only unique up to a shift by a constant. In particular, the condition on $T'(h)$ cannot be used to determine this constant of integration. Once the constant is chosen to satisfy the condition on $T'(0)$, the condition on $T'(h)$ is automatically satisfied by the "compatibility condition" from the previous paragraph.