DISCLAIMER: This reply only handles integration by parts to show how to obtain the Neuman condition for this problem. To obain the axi-symmetric variational formulation, a change of integration coordinate system has to be applied as well, see the old Q&A-forum
So, you have the equation:
C\partial T/\partial t - \partial (K \partial T)/\partial r = q .
Using finite difference temporal discretization and Integration by parts yields
(C (T^{i+1}- T^i)/\delta t, v)_\Omega + K(\partial T^{i +1}/\partial r, \partial v/\partial r)_\Omega - (K n\partial T^{i +1}/\partial r, v )_{\partial\Omega}=(q,v)_\Omega, where n is ±1 depending on if you are at x=0 or x=1.
Since you have a zero neumann condition, the boundary terms falls away and you are left with volume terms.
To get the point source q, you can follow Help a starter! on how to make a point source (in that thread it has the name f).