Heat Equation part 2 a slight modification

While writing the code for the previous post I slightly modified the code in order to add 2 ‘peaks of heat’. Now, by looking at what I wrote I am not quite sure that the solution is consistent (physically speaking). At first glance we can see that there are three ‘heat spots’. If you imagine of heating the rod in two spots and cooling it in between, the situation is the following:



By setting L = 3*pi and making some other small changes we obtain:

Physically we could interpret the result as follows:
the heat flows spontaneously towards areas where the temperature is lower, therefore the temperature should go up in the middle of the rod (since it has been cooled down and has a lower temperature compared its surroundings) and go down in the sides. Assuming that the left and right boundaries of the rod are kept at constant temperature through an ideal thermostat.

	import numpy as np
	from numpy import pi
	import matplotlib.pyplot as plt
	import matplotlib.animation as animation
	fig = plt.figure()
	fig.set_dpi(100)
	ax1 = fig.add_subplot(1,1,1)
	#Diffusion constant
	k = 2
	#Scaling factor
	scale = 5
	#Length of the rod (0,L) on the x axis
	L = 3*pi
	#Initial contitions
	x0 = np.linspace(0,L,10000)
	t0 = 0
	temp0 = 5 #Temperature of the thermostat
	#Increment
	dt = 0.02
	#Heat equation
	def u(x,t):
	return temp0 + scale*np.exp(-k*t)*np.sin(x)
	#Gradient of the Heat Equation
	def grad_u(x,t):
	#df/dx #df/dt
	return scale*np.array([np.exp(-k*t)*np.cos(x),-k*np.exp(-k*t)*np.sin(x)])
	a = []
	for i in range(500):
	value = u(x0,t0) + grad_u(x0,t0)[1]*dt
	t0 = t0 + dt
	a.append(value)
	k = 0
	def animate(i):
	global k
	x = a[k]
	k += 1
	ax1.clear()
	plt.plot(x0,x,color='red',label='Temperature while cooling after heating')
	plt.plot(x0,temp0*np.ones(x0.shape),color='green',label='Thermostat temperature of the rod')
	plt.grid(True)
	plt.legend()
	plt.ylim([-4,2.5*scale])
	plt.xlim([-0.2,L+0.2])
	anim = animation.FuncAnimation(fig,animate,frames=360,interval=20)
	plt.show()

view raw heatEquation2.py hosted with ❤ by GitHub

Why not  checking if this solution is mathematically and physically consistent?