You are correct in that you use Neumann BC at the edge to express thermal insulation. Because you work in polars, this means that
$$\left[ \frac{\partial T}{\partial r}\right]_{r=b} = 0$$
Keep in mind that Laplaces's equation is for the steady-state temperature. For a changing temmperature, you would be solving the heat equation, which has a time derivative in it.
For steady state temp, you should be able to show that
$$T(r,\theta) = \left(\frac{A}{r}+ B r \right )\cos{\theta}$$
There are no higher-order terms because of the nature of your boundary conditions. Now $T(a,\theta) = \cos{\theta}$ so that
$$\frac{A}{a} + B a = 1$$
Your second boundary condition on thermal insulation implies that
$$-\frac{A}{b^2} + B = 0$$
For (iii), the heat flux is just the derivative with respect to $r$ at the inner surface, multiplied by $-\kappa$.
For the steady state, you are solving
$$u_t=0 \implies u_{xx}=-f(x)$$
The general solution to this, given $f$, is pretty straightforward:
$$u(x) = \begin{cases} c_1 x+c_2 & x \in \left (0,\frac{L}{2}\right)\\-\frac12 H x^2 + d_1 x+d_2 & x \in \left (\frac{L}{2},L\right) \end{cases} $$
Yes, we have different constants in each branch of the integration interval. We thus need four conditions. Two come from the boundary conditions:
$$u(0)=u(L)=0$$
The other two come from continuity requirements for $u$ and its derivative (otherwise, the heat flow is discontinuous, not what should happen in a steady-state solution). Thus,
$$c_1 \frac{L}{2}+c_2 = -\frac18 H L^2 + d_1 \frac{L}{2} + d_2$$
$$c_1 = -\frac12 H L + d_1$$
Four equations, four unknowns, solve. It's not as bad as you may think. For example, $c_2=0$. The rest fall out from simple substitutions. The result I get is
$$u(x) = \begin{cases} \frac18 H L x & x \in \left (0,\frac{L}{2}\right)\\-\frac12 H x^2 + \frac{5}{8} H L x-\frac18 H L^2 & x \in \left (\frac{L}{2},L\right) \end{cases} $$
EDIT
To see what the above temperature looks like, rewrite the above result in the following form:
$$f(y) = y \,\theta \left ( \frac12-y\right) - (4 y-1)(y-1)\, \theta \left (y- \frac12\right) $$
where $\theta$ is the Heaviside step function, $y=x/L$, and $f(y) = 8 u(L y)/(H L^2)$.
Here is a plot of $f(y)$:
Best Answer
For simplicity, let's suppose that your disk is the unit disk. The heat equation is $u_t = k\Delta u$. Steady state means that the temperature $u$ does not change; thus $u_t=0$ and you are left with Laplace's equation: $\Delta u=0$ subject to $u(1,\theta)=f(\theta)$. The solution may then be written:
$$u(r,\theta )=\frac{a_0}{2}+\sum _{n=1}^{\infty } r^n\left(a_n\cos (n \theta )+b_n\sin (n \theta )\right),$$
where $a_n$ and $b_n$ are the Fourier coefficients of $f$, as justified below. Alternatively, the $a_n$s and $b_n$s may be replaced with their integral representations, the order of summation and integration may be flipped, and (after a lot of simplification) this morphs into an integral representation the Poisson integral formula:
$$u(r,\theta )=\frac{1}{2\pi }\int _{-\pi }^{\pi }f(\phi )\frac{R^2-r^2}{R^2+r^2-2R r \cos (\theta -\phi )}d\phi.$$
This is particularly nice for numerical computation.
An example
As a cute example, suppose that $f(\theta )=\sin (2\theta )$. Then, $f$ is already a very short Fourier series and the solution is $u(r,\theta )=r^2\sin (2\theta )$. This looks something like so:
Note that $\sin(2\theta)$ meanders up and down on the boundary and the solution obeys many of the properties that we expect from a steady state solution. In particular, the Maximum Principle is satisfied. This is probably how the Pringle was discovered.
Separation of variables
Here is the outline of the separation of variables required to derive the Fourier series representation. We start with the polar specification of the problem.
$$u_{rr}+\frac{1}{r}u_r+\frac{1}{r^2}u_{\theta \theta }=0,\text{ }u(1,\theta )=f(\theta ).$$
Setting $u(r,\theta )=R(r)\Theta (\theta )$, the PDE becomes
$$R^{\prime\prime }\Theta +\frac{1}{r}R'\Theta +\frac{1}{r^2}R\, \Theta ^{\prime\prime }=0$$
which separates into
$$-\frac{r^2R^{\prime\prime }+r\, R'}{R}=\frac{\Theta ^{\prime\prime }}{\Theta }=-\lambda$$
leading to the two ODEs
$$r^2R^{\prime\prime }+r\, R'=\lambda \, R\text{ }\text{and}\text{ }\Theta ^{\prime\prime }=-\lambda \, \Theta .$$
To find the eigenstructure, we focus on the $\theta$ equation. First, note that $\lambda =0$ is an eigenvalue with any constant a representative eigenfunction. Now, implicit in the $\theta$ equation are periodic boundary conditions $\Theta (0)=\Theta (2\pi )$ and $\Theta '(0)=\Theta '(2\pi )$. Furthermore, the general solution of the $\theta$ equation is
$$\Theta (\theta )=a \cos \left(\sqrt{\lambda }\theta \right) + b \sin \left(\sqrt{\lambda }\theta \right)$$
which has derivative
$$\Theta '(\theta )=-a\sqrt{\lambda } \sin \left(\sqrt{\lambda }\theta \right) + b\sqrt{\lambda } \cos \left(\sqrt{\lambda }\theta \right).$$
Thus, the boundary conditions become
$$a =a \cos \left(2\pi \sqrt{\lambda }\right) + b \sin \left(2\pi \sqrt{\lambda }\right)$$
and
$$b\sqrt{\lambda }=-a\sqrt{\lambda } \sin \left(2\pi \sqrt{\lambda }\right) + b\sqrt{\lambda } \cos \left(2\pi \sqrt{\lambda }\right).$$
This is equivalent to the system
$$\left(\cos \left(2\pi \sqrt{\lambda }\right)-1\right)a + \sin \left(2\pi \sqrt{\lambda }\right)b = 0\\ \\ a \sin \left(2\pi \sqrt{\lambda }\right) + \left(1-\cos \left(2\pi \sqrt{\lambda }\right)\right)b = 0$$
or
$$\left( \begin{array}{cc} \left(\cos \left(2\pi \sqrt{\lambda }\right)-1\right) & \sin \left(2\pi \sqrt{\lambda }\right) \\ \sin \left(2\pi \sqrt{\lambda }\right) & \left(1-\cos \left(2\pi \sqrt{\lambda }\right)\right) \\ \end{array} \right)\left( \begin{array}{c} a \\ b \\ \end{array} \right)=\left( \begin{array}{c} 0 \\ 0 \\ \end{array} \right).$$
This has a non-trivial solution if and only if the matrix has determinant zero, i.e.
$$\left(\cos \left(2\pi \sqrt{\lambda }\right)-1\right)\left(1-\cos \left(2\pi \sqrt{\lambda }\right)\right) - \sin ^2\left(2\pi \sqrt{\lambda }\right) = 0.$$
Using a few trig identities, this simplifies down to $-4\sin ^2\left(2\pi \sqrt{\lambda }\right)=0$. This implies that $\sqrt{\lambda}$ must be a positive integer or $\lambda =n^2$ for some $n\in \mathbb{N}$.
Next, we must solve the $r$ equation. Setting $\lambda =n^2$ in the $r$ equation we get
$$r^2R^{\prime\prime }+r\, R'=n^2R.$$
This is called a Cauchy-Euler equation and has general solution $R(r)=c_nr^{-n}+d_nr^n$, which is quite easy to check. Since the solution should be bounded, we must have $c_n=0$ for all $n$. As a result, we have $R_n(r)=d_nr^n$ and
$$u(r,\theta )=\frac{a_0}{2}+\sum _{n=1}^{\infty } r^n\left(a_n\cos (n \theta )+b_n\sin (n \theta )\right).$$
Note that the $a_0$ arises from the zero eigenvalue and that the $d_n$ can be sucked up into the other constants.
Finally, we need to choose the $a_n$s and $b_n$s so that the boundary condition is satisfied, i.e. we want
$$u(1,\theta )=\frac{a_0}{2}+\sum _{n=1}^{\infty } \left(a_n\cos (n \theta )+b_n\sin (n \theta )\right) = f(\theta).$$
Thus, the $a_n$s and $b_n$s are just the full Fourier coefficients of $f$ as claimed.