I am trying to solve the heat equation for a system made of a cylinder (hot) that is suddenly immersed in a cooling medium (cold). I am doing this using a 1D approximation with single differential equation (in cylindrical coordinates) and physical properties that vary in space. The code is:

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