diamcyl = 0.800; T1 = 140; T2 = -20;

k2 = 0.58; (* W/mK ) rho2 = 1000; ( kg/m3 ) Cp2 = 4200; ( J/KgK *)

k1 = 0.128; (W/(mK)) rho1 = 800; (* kg/m3 *) Cp1 = 1670; (J/(KgK))

(properties as a function of space) rho[x_] := (rho1 + (rho2 - rho1)* UnitStep[x - diamcyl/2]); k[x_] := 10^6*(k1 + (k2 - k1)* UnitStep[x - diamcyl/2]);
Cp[x_] := (Cp1 + (Cp2 - Cp1)* UnitStep[x - diamcyl/2]);

heateq = 1/xD[xk[x]*D[u[x, t], x], x] == rho[x]*Cp[x]*D[u[x, t], t];

u0[x_] := (T1 + (T2 - T1)*UnitStep[x - diamcyl/2]); solm20 = First[ NDSolve[{heateq, u[x, 0] == u0[x], u[10, t] == T2, (D[u[x, t], x] /. x -> 0.001) == 0}, u, {x, 0.001, 2}, {t, 0, 5}, Method -> {"MethodOfLines", "SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> 1000}}]];

The solution transient temperature profile that I obtain looks smooth and is continuous over the boundary of the cylinder. However, the heat flux exhibits a discontinuity over this boundary.

Show[Table[ Plot[u[x, t] /. solm20, {x, 0.001, 2}, PlotRange -> {{0, 1}, {T2, T1}}, PlotRangePadding -> {None, 10}, Prolog -> {LightGray, Rectangle[{0, -100}, {diamcyl/2, 200}]}, GridLines -> {{diamcyl/2}, {T2}}, FrameLabel -> {"distance [mm]", "temperature [[Degree]C]"}], {t, 0, 2, 0.05}],

Plot[u0[x], {x, 0.01, 2}, PlotRange -> All, Exclusions -> None, PlotStyle -> Red]]

Show[Table[ Plot[Evaluate[k[x]*D[u[x, t], x] /. solm20], {x, 0.01, 2}, PlotRange -> {{0, 1}, All}, PlotRangePadding -> {None, 10}, GridLines -> {{diamcyl/2}, {T2}}, FrameLabel -> {"distance [mm]", "temperature [[Degree]C]"}], {t, 0.01, 2, 0.05}]]

How to obtain a solution with continuous heat flux? Or, I could also ask it in this way: how to impose a Neumann condition over the boundary of the cylinder? Thanks

diamcyl = 0.800;
T1 = 140;
T2 = -20;

k2 = 0.58;    (* W/mK *)
rho2 = 1000;   (* kg/m3 *)
Cp2 = 4200;  (* J/KgK *)

k1 = 0.128; (*W/(mK)*)
rho1 = 800; (* kg/m3 *)
Cp1 = 1670; (*J/(KgK)*)

(*properties as a function of space*)
rho[x_] := (rho1 + (rho2 - rho1)* UnitStep[x - diamcyl/2]);
k[x_] := 10^6*(k1 + (k2 - k1)* UnitStep[x - diamcyl/2]);    
Cp[x_] := (Cp1 + (Cp2 - Cp1)* UnitStep[x - diamcyl/2]);

heateq = 1/x*D[x*k[x]*D[u[x, t], x], x] == rho[x]*Cp[x]*D[u[x, t], t];

u0[x_] := (T1 + (T2 - T1)*UnitStep[x - diamcyl/2]);
solm20 = First[
   NDSolve[{heateq, u[x, 0] == u0[x], 
     u[10, t] == T2, (D[u[x, t], x] /. x -> 0.001) == 0}, 
    u, {x, 0.001, 2}, {t, 0, 5},
    Method -> {"MethodOfLines", 
      "SpatialDiscretization" -> {"TensorProductGrid", 
        "MaxPoints" -> 1000}}]];

The solution transient temperature profile that I obtain looks smooth and is continuous over the boundary of the cylinder. However, the heat flux exhibits a discontinuity over this boundary.

Show[Table[
  Plot[u[x, t] /. solm20, {x, 0.001, 2}, 
   PlotRange -> {{0, 1}, {T2, T1}}, PlotRangePadding -> {None, 10}, 
   Prolog -> {LightGray, Rectangle[{0, -100}, {diamcyl/2, 200}]}, 
   GridLines -> {{diamcyl/2}, {T2}}, 
   FrameLabel -> {"distance [mm]", "temperature [\[Degree]C]"}], {t, 
   0, 2, 0.05}],
 
 Plot[u0[x], {x, 0.01, 2}, PlotRange -> All, Exclusions -> None, 
  PlotStyle -> Red]]

Show[Table[
  Plot[Evaluate[k[x]*D[u[x, t], x] /. solm20], {x, 0.01, 2}, 
   PlotRange -> {{0, 1}, All}, PlotRangePadding -> {None, 10}, 
   GridLines -> {{diamcyl/2}, {T2}}, 
   FrameLabel -> {"distance [mm]", "temperature [\[Degree]C]"}], {t, 
   0.01, 2, 0.05}]]

How to obtain a solution with continuous heat flux? Or, I could also ask it in this way: how to impose a Neumann condition over the boundary of the cylinder?