I’m attempting to model a situation in which a polymer initially at a higher temperature is sandwiched between two cooler metallic mold pieces with conductive boundary conditions in between the polymer and mold. On either outer side of the mold, there is convection with another fluid. Since the boundary conditions at either side of the sandwiched polymer where contact is made with the mold are the same, I’m focusing on only one of the mold pieces, with the mold spanning from $0 leq z leq d$ and the polymer spanning from $-w leq z leq 0$. The following equations model the setup:

$$frac{partial T_{mold}(z,t)}{partial t}=a_{mold}frac{partial^2T_{mold}(z,t)}{partial z^2}, 0 leq z leq d$$

$$frac{partial T_{poly}(z,t)}{partial t}=a_{poly}frac{partial^2T_{poly}(z,t)}{partial z^2}, -w leq z leq 0$$

Boundary conditions:

$$k_{mold}frac{partial T_{mold}(z,t)}{partial z}|_{z=d}=h_c(T_{mold}(z,t)-T_c)$$

$$-k_{mold}frac{partial T_{mold}(z,t)}{partial z}|_{z=0}=-k_{poly}frac{partial T_{poly}(z,t)}{partial z}|_{z=0}$$

$$-k_{mold}frac{partial T_{mold}(z,t)}{partial z}|_{z=-w}=-k_{poly}frac{partial T_{poly}(z,t)}{partial z}|_{z=-w}$$

Initial conditions:

$$T_{mold}(z,t=0)=T_m$$

$$T_{poly}(z,t=0)=T_p$$

I tried to approach solving this based on the wonderful iterative approach of @Alex Trounev from this post using the following code (the signs on NeumannValue might be incorrect due to my inexperience):

```
amold = 0.004;
Tm = 323;
apoly = 0.0001;
Tp = 433;
Tc = 303;
h = 5000;
d = 0.005;
w = 0.005;
poly(0)(z_, t_) := Tp;
Do(mold =
NDSolveValue({D(Tmold(z, t), t) -
amold*D(Tmold(z, t), {z, 2}) ==
NeumannValue(-h (Tmold(z, t) - Tc)/
kmold, z == d) +
NeumannValue(kpoly/kmold)*(D(poly(i - 1)(z, t), z) /. z -> 0), z == 0),
Tmold(z, 0) == Tm},
Tmold, {z, 0, d}, {t, 0, 600},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"FiniteElement",
"MeshOptions" -> {"MaxCellMeasure" -> {"Length" -> 0.001}}}});
poly(i) =
NDSolveValue({D(Tpoly(z, t), t) -
apoly*D(Tpoly(z, t), {z, 2}) ==
NeumannValue(kmold/kpoly)*(D(mold(z, t), z) /. z -> 0), z == 0) +
NeumannValue(kmold/kpoly)*(D(mold(z, t), z) /. z -> 0), z == -w),
Tpoly(z, 0) == Tp},
Tpoly, {z, -w, 0}, {t, 0, 600},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"FiniteElement",
"MeshOptions" -> {"MaxCellMeasure" -> {"Length" ->
0.001}}}});, {i, 1, 20})
```

However, I get a very unstable solution, and the kernel crashes after ~30 iterations. I also tried solving the system of PDEs without an iterative approach below, but I get the error that the boundary condition at $z=-w$ is not specified on a single edge of the boundary of the computational domain:

```
sys1 = {D(Tmold(z, t), t) -
amold*D(Tmold(z, t), {z, 2}) ==
NeumannValue(-h (Tmold(z, t) - Tc)/
kmold, z == d),
Tmold(z, 0) == Tm,
D(Tpoly(z, t), t) ==
apoly*
D(Tpoly(z, t), {z, 2}), -kmold*(D(Tmold(z, t), z) /. z -> 0) == -kpoly*(D(Tpoly(z, t), z) /. z -> 0), -kmold*(D(Tmold(z, t), z) /. z -> 0) == -kpoly*(D(Tpoly(z, t), z) /. z -> -w),
Tpoly(z, 0) == Tp};
heat = NDSolveValue(sys1, {Tmold,Tpoly}, {z, -w, d}, {t, 0, 600},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"FiniteElement",
"MeshOptions" -> {"MaxCellMeasure" -> {"Length" -> 0.001}}}})
```

I was wondering if anyone had insight into how I can improve either of these two approaches. Thank you again in advance for any help!