a = 0.5*10^-3; b = 20*10^-3; \[Rho] = 11.3*10^3; c = 130; \[Lambda] = 300; i = 10; R = 0.01; H = 15; fun = \[Rho]*c*D[T[r, z, t], t] - \[Lambda]* Laplacian[T[r, z, t], {r, z}] - (i^2 R)/(\[Pi]*a^2*b); ic = T[r, z, 0] == 300.15; bc1 = fun == NeumannValue[H*(T[r, z, t] - T[r, z, 0]), z == 0]; bc2 = fun == NeumannValue[-H*(T[r, z, t] - T[r, z, 0]), r == a]; bc3 = fun == NeumannValue[-H*(T[r, z, t] - T[r, z, 0]), z == b]; sol = NDSolveValue[{bc2, bc1, bc3, ic}, T[r, z, t], {r, 0, a}, {z, 0, b}] Clear["Global`*"]