Consider a hollow cylinder of different outer radius and inner radius and two different temperatures are maintained at the outer and inner surfaces such that inner temperature is higher. Because of this heat will flow radially outward.
Now I don't understand why we need to use integration to find rate of heat flow that is, why isn't the temperature gradient constant? Why can't we directly apply heat conduction formula?
Okay, the "cylinder" thing means that we can think of this as a 2D system; we can assume that the temperature of the cylinder has reached a steady state. Then there are three constraints at play here:
You are right that if you have two boxes separated by a sheet, then the linear transfer equilibrium is solved for a linear temperature gradient. That is effectively a 1D problem.
But when you generalize this to 2D and higher dimensions, the linear temperature gradient doesn't give you this linear transfer equilibrium. Why not? Well think about the absolute simplest case, two boxes at T=0 and T=1 which are at opposite ends of a 3x3 matrix separated by walls. It's not too hard to work out that the temperature solution of this:
+-----+-----+-----+
| T=0 | | |
+-----+-----+-----+
| | | |
+-----+-----+-----+
| | | T=1 |
+-----+-----+-----+
does not look linear:
+-----+-----+-----+ /|\T *
| 0 | 1/3 | 1/2 | |
+-----+-----+-----+ | *
| 1/3 | 1/2 | 2/3 | | *
+-----+-----+-----+ | *
| 1/2 | 2/3 | 1 | | distance along
+-----+-----+-----+ +-*--------------> diagonal
So the way to think of this example is that, say, focusing on the 1/3 cell, you see a flow of 1/3 "out" to the 0 cell exactly balanced by a flow of 2 * 1/6 "in" from the 1/2 cells next to it. Notice that we're missing two values, 1/6 and 5/6, needed for a "linear progression" across the diagonals. We instead see a graph which has a little bit of a curve when we look along the diagonal. There's an even-simpler case where we force the 3 cells in the top corner to be 0, and that's solved by:
+-----+-----+-----+
| 0 | 0 | 1/4 | /|\
+-----+-----+-----+ | *
| 0 | 1/4 | 1/2 | |
+-----+-----+-----+ | *
| 1/4 | 1/2 | 1 | | *
+-----+-----+-----+ +-*--*----------->
Again, the progression is not linear! {0, 1/4, 1/2, 1} is missing a value 3/4 which is nowhere to be found.
The bottom line is: surface area matters. When you've got this increasing/decreasing surface area, even in the "transfer region" you don't necessarily see linear gradients, because it's not just how much energy is in the box that matters -- it's how many boxes are nearby, that this one is flowing into.
OKAY. Returning to your case, how do we solve the conjunction of (2) and (3) for a circular symmetry? The answer is actually really simple. First we notice that the "local" heat flux at any point $J$ must be dropping like $J = k / r$ because, by (1), it is constant along the circles we're drawing, so the integral of it is simply $2 \pi r J$. If that's constant with respect to radius then we have $J \propto 1/r$.
Now we use (3). The gradient in the radial direction is actually still just d/dr, so we just have to integrate $k/r$ to get: $T(r) ~=~ A ~-~ B~ \operatorname{ln} (k / r) $. Those look like three independent constants but because of log-laws we know that they are not, so we can just set $k = R_0$ or so. Solve for A and B such that $T(R_0) = T_0$ and $T(R_1) = R_1$, if you like.
Finally notice that your intuition is right if the cylinder gets super big and thin, in which case we can write $r = R + \delta r$ and expand to first order, getting a linear temperature gradient.
The fact that it is a cylinder has to be taken into account. Usually $\frac{dQ}{dt} = \lambda \frac{A}{d}\Delta T$, where A is the surface perpendicular to the heat flow and d is the thickness of the material between the to temperatures. As the temperatures are being maintained, the heatflow is a constant, but being a cylinder the outer surface is larger than the inner surface. So you have to work with polar coordinates. This leads to $$-\triangledown^2 T=\frac{\delta^2 T}{\delta r^2}-\frac{1}{r}\frac{\delta T}{\delta r} $$ So the integrations stems from the fact that it is a cylinder.
You will find more information in this article: INTRODUCTION TO ENGINEERING HEAT TRANSFER
and here: Heat Transfer/Conduction
The temperature does not depend linearly on the radius because the flow lines of heat flux are not parallel. This means that the heat flux, and therefore the temperature gradient, decrease as the flux lines diverge, and increase as the converge. Here the flow lines are radial.
By conservation of energy, the rate of heat passing through any cylindrical surface concentric with the inner and outer surface (i.e. the surface is everywhere normal to the heat flow lines) must be the same for all such surfaces at thermal equilibrium (i.e. steady state not to be confused with thermodynamic equilibrium), otherwise, heat would be gathering or "pooling" between two cylindrical surfaces with different heat fluxes through them, to fulfill conservation of energy.
But the area of such a cylindrical surface is proportional to its radius. Therefore, the heat flux per unit area must decrease in inverse proportion to the cylinder's radius to ensure that the total heat flux through each such cylindrical surface (in joules per second) is constant. But the heat flux is proportional to the the temperature gradient $\mathrm{d}T/\mathrm{d}r$. Therefore, we have:
$$r\,\frac{\mathrm{d}T}{\mathrm{d}r} = const$$
or
$$T = A\,\log r + B$$
and you match this to your boundary conditions $T(r_i) = T_i,\,T(r_o)=T_o$ with the appropriate choice of $A$ and $B$.
