NDSolve solution

Clear[f];
f[r_]  =  Which[r < 1, 0, r < 2, 1, True, 0];
\(\(g1 = 
        ParametricPlot3D[
	{r\ Cos[phi], r\ Sin[phi], f[r]}, 
	{phi, 0, 
            2\ Pi}, {r, 0, 3}];\)\)

Remove[T]
\(\(nsol = 
        NDSolve[
	{D[r\ D[T[r], r], r] == \(-\ r\)\ *f[r], 
		T[3] == 
              0, \(\(Derivative[1]\)[T]\)[0.01] == 0}, 

          T[r], {r, 0.01, 3}, 
	Method\  -> \ RungeKutta];\)\)
Tn[r_]  =  T[r] /. nsol[[1]];

\(\(\(g2 = Plot[Tn[r], {r, 0.01, 3}, PlotRange\  -> \ All];\)\(
\)
    \( (*CreateEPS["heatdistribution", g2]*) \)\)\)

\(\(ParametricPlot3D[
	{r\ Cos[phi], r\ Sin[phi], Tn[r]}, 
	{phi, 0, 
          2\ Pi}, {r, 0.01, 3}];\)\)

Clear[T, Tn]
\(\(nsol = 
        NDSolve[
	{D[r\ D[T[r], r], r] == \ r\ *f[r], 
		T[3] == 
              0, \(\(Derivative[1]\)[T]\)[0.01] == 0}, 

          T[r], {r, 0.01, 3}];\)\)
Tn[r_]  =  T[r] /. nsol[[1]];
Plot[Tn[r], {r, 0.01, 3}, PlotRange  ->  All];

Exact Solution

Clear[fe, dfe, c1, c2, c3, c4, x]

\(\(fe[x_]\  := \ 
        Which[x < 1, \ c1, x < 2, \(-x^2\)/4 + c3\ Log[x] + c4, 
		x < 4, 
          c2\ Log[x/3]];\)\)
dfe[x_]  =  D[fe[x], x];

{fe[0.99], fe[1]}
{dfe[0.99], dfe[1]}

{fe[1.99999], fe[2]}
{dfe[1.99999], dfe[2]}

\({c1, c2, c3, c4} = 
	
      LinearSolve[{{1, 0, 
              0, \(-1\)}, 
			\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {0, 
              0, \(-1\), 0}, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 	{0, 
              Log[2/3], \(-Log[
                  2]\), \(-1\)}, 
			\ \ \ \ \ \ \ \ \ \ \ \ \ \ {0, 
              1/2, \(-1\)/2, 0}}, 
{\(-1\)/4, \(-1\)/2, \(-1\), \(-1\)}] // 
        N\)

\(\(fen[x_]\  := \ 
        Which[x < 1, \ c1, x < 2, \(-x^2\)/4 + c3\ Log[x] + c4, 
		x < 4, 
          c2\ Log[x/3]];\)\)
gexact = Plot[fen[x], {x, 0, 3}];

FEM Solution

<< BVP.m;

Clear[a, b, c, f, x, t, n];
a = Function[x, x];
b = Function[x, 0];
f = Function[x, Which[x < 1, 0, x < 2, -x, True, 0]];
n = 3;
x = Table[t, {t, 0, 3, 3/n}];
data = BVP[a, b, f, x, {"N", "D"}, {0, 0}];
\(\(g1 = 
        ListPlot[
          data, 
		PlotStyle\  -> \ 
            PointSize[0.02], 
		DisplayFunction\  -> \ 
            Identity];\)
\)
g2 = Show[gexact, g1];

<< Interpol2.m

fn[x_]  :=  Interpol2[x, data]
Plot[{fn[x] - fe[x]}, {x, 0, 3}, PlotRange -> All];

\(res = 
      Table[{n = 3*k, x = Table[t, {t, 0, 3, 3/n}]; 

          data = BVP[a, b, f, x, {"N", "D"}, {0, 0}]; 
          data[\([1, 2]\)] - fe[0]}, {k, 1, 30, 3}]\)

ListPlot[res, PlotJoined  -> True, PlotRange -> All];

\(\(ListPlot[Log[res], 
	PlotJoined\  -> True, 

        PlotRange -> All];\)\)

lres = Log[res];
n = Length[res];
(lres[[1, 2]] - lres[[n, 2]])/(lres[[1, 1]] - 
          lres[[n, 1]])