Does anybody know how I should deal with this type of problem?
Transform the PDE in homogenous nondimensional form and plot steady-state profiles for both cases.
Should I find a transformation to convert the above PDE to homogeneous ones? is there any transformation to convert non-homogeneous (like the above) to the homogenous?

I think the question does not ask for a transformation that removes the source term $q$. Maybe the word 'homogeneous' refers to homogeneity of dimensions and units; a property that should be verified by the dimensionless form.
To do so we introduce dimensionless quantities $$ \tau^* = \frac{\tau}{\tau_c}, \quad z^* = \frac{z}{z_c}, \quad T^* = \frac{T}{T_c}, \quad q^* = \frac{q}{q_c} , $$ where the quantity with index $c$ is a characteristic quantity (or reference quantity). The differential operators satisfy $$ \frac{\partial}{\partial \tau} = \frac{\partial \tau^*}{\partial \tau}\frac{\partial}{\partial \tau^*} = \frac{1}{\tau_c}\frac{\partial}{\partial \tau^*} , $$ etc. Thus the PDE rewrites as $$ \frac{\partial T^*}{\partial \tau^*} = D \frac{\tau_c}{{z_c}^2}\frac{\partial^2 T^*}{\partial {z^*}^2} + \frac{D}{k} \frac{q_c\tau_c}{T_c} q^* . $$ where $D = {k}/({\rho C_p})$ is the diffusivity constant. Now, we look for appropriate choices of the characteristic quantities that remove the parameters of the equation (all coefficients equal one). Therefore, we can write two equations of four unknowns $\tau_c$, $z_c$, $T_c$ and $q_c$. Here, let us impose the characteristic length $z_c = 1$ m of the spatial domain and the characteristic temperature $T_c = 1$ K. Thus, we have $$ \tau_c = \frac{{z_c}^2}{D}, \quad q_c = \frac{kT_c}{{z_c}^2} , $$ and the PDE rewrites as $$ \frac{\partial T^*}{\partial \tau^*} = \frac{\partial^2 T^*}{\partial {z^*}^2} + q^* . $$ Now, it remains to compute the steady-state solution by solving the differential equation $T^*_{z^*z^*} = -q^*$ with the prescribed dimensionless boundary conditions and source terms.