Hi! I'm a beginner level user of Mathematica, and still learning it. Basically, I'm trying to monitor the natural cooling of an imaginary plastic slab with the top/bottom side insulated, only the left/right side is dissipating heat to the environment. Thus the model reduced to a 1D heat transfer model. The slab has a thickness of 10mm, and initial uniform temperature of 275C, the environment temperature is 25C. I'm trying to implement a Neumann type boundary condition at both left (x = -5) and right (x = 5) side, i.e., k*dT/dx = h(T-Tenvironment) at boundary, but I keep getting error messages that the boundary and initial conditions are inconsistent and more importantly the results looks unstable, my code looks like this, can anyone give me a little help? Would be really appreciated!
h = 13/(1000^2);(*units are converted to mm*)
k = 0.3/1000;
c = 1030;
rou = 1380/(1000^3);
heq = rou*c*D[u[t, x], t] == k*D[u[t, x], {x, 2}];
ic = u[0, x] == 275;
bcs1 = (D[u[t, x], x] /. x -> 5) + h*(u[t, 5] - 25) == 0;
bcs2 = (D[u[t, x], x] /. x -> -5) + h*(u[t, -5] - 25) == 0;
sol = NDSolve[{heq, ic, bcs1, bcs2}, u[t, x], {t, 0, 20}, {x, -5, 5}];
F[t_, x_] = u[t, x] /. sol[[1]];
Plot3D[F[t, x], {x, -5, 5}, {t, 0, 15}, PlotRange -> All,
AxesLabel -> Automatic]
Attachments: