I’m trying to solve following PDE (heat equation):

$ $ \begin{cases} u_t = a \, u_{xx} \ u(x,0) = 0\ \lim_{x\to \infty}u(x,t) =0\ \alpha\, (\theta_0-u(0,t))+\dot{q}_0=-\lambda u_x(0,t) \end{cases}$ $

Where basically I have an initial temperature of $ 0$ everywhere, a constant heat flux at the beginning of the rod, and convection between air and the rod at its beginning ($ \theta_0$ is the air temperature which is assumed to be constant).

I found following analytical solution for the problem:

$ $ u(x,t) = \frac{\dot{q}_0+\alpha \, \theta_0}{\alpha}\left[ \mathrm{erfc}\, \left(\frac{x}{2\sqrt{a \, t}} \right) -\mathrm{exp}\,\left(\frac{\alpha}{\lambda}x+a \frac{\alpha^2}{\lambda^2}t \right)\mathrm{erfc}\,\left(\frac{x}{2\sqrt{a\, t}} + \frac{\alpha}{\lambda} \sqrt{a\, t} \right) \right] $ $

which is physically meaningful. With mathematica, however, I get some meaningless results, probably due to my not that good mathematica skills.

This is what I tried to do:

`λ = 0.8; c = 880; ρ = 1950; a = λ/(c ρ) α = 15; θair = 0; heqn = D[u[x, t], t] == a D[u[x, t], {x, 2}]; ic1 = u[x, 0] == 0; ic2 = α (θair - u[0, t]) + 650 == -λ Derivative[1, 0][u][0, t]; sol = DSolveValue[{heqn, ic1, bc1}, u[x, t], {x, t}][[1, 1, 1]] `

Which leads me to a complex solution (plotted below)

`DensityPlot[sol, {t, 0, 2*3600}, {x, 0, 0.1}, PlotLegends -> Automatic, FrameLabel -> {t[s], x[m]}] `

[

`Plot[Evaluate[Table[sol, {t, 3600, 7200, 3600}]], {x, 0, .1}, PlotRange -> All, Filling -> Axis, Axes -> True, AxesLabel -> {x[m], θ[°C]}] `

I’m aware of the fact that I didn’t consider the condition at infinity. To do this, I tried to follow this answer without success. Also, mathematica finds the solution without this condition as soon as $ \alpha = 0$ .

This is the plot of the analytical solution I get:

`Plot[{u[x, 600], u[x, 3600], u[x, 7200]}, {x, 0, .2}, Filling -> Axis, Axes -> True, AxesLabel -> {x[m], θ[°C]}] `