Spike in Temperature and Dielectric Constant Loop Plot

Hi I was wondering if I can get some advice as to why there is a large spike in the plot of my code below. I am trying to plot the temperature and dielectric constant of a 1D heated slab (heated on the left side) as a function of x. The dielectric constant is related to T with 2 different equations above and below the critical temperature TC.

I suspect the spike may be due to very small values of TNew-TC in some places, but the value of the spike seems to increase with more iterations which doesn't make sense. I was wondering if I was missing something or if there is a bug I have missed?

Thanks in advance.

%Initialise Variables:

tNew=120 %seconds, simulation time

dt=1 %seconds, timestep

L=0.05; %metres, slab thickness

nx=20; %metres, number of steps in x direction

dx=L/nx;%metres, size of step in x direction

x=dx/2:dx:L-dx/2;%metres, position in x direction

k=6 %thermal conductivity

rho=6020 %density

C=527 %specific heat capacity

const=(k)/(rho*C) %thermal diffusivity

T0=298; %Kelvin, initial temperature of slab

TW=450; %Kelvin, initial tempature of left surface

TE=298; %Kelvin, initial tempature of right surface

TC=ones(nx,1)*408.15; %Kelvin, first order transition temperature

K=1.6*(10^5)%Curie constant

%Create Matrix with Initial values:

T=ones(nx,1)*T0;

TNew=zeros(nx,1);

v=VideoWriter('myVideo')

v.FrameRate=5

open(v)

t=0:dt:tNew

time = -1

for i = 1:length(t)

for j = 2:nx-1

TNew(j) = T(j) + dt*const*((T(j+1)-T(j))/dx^2 + (T(j-1)-T(j))/dx^2);% Master equation

if (TNew(j)>=TC(j))

ENew(j) = K./(TNew(j)-TC(j))

else

ENew(j) = K/(2*(TC(j)-TNew(j)))

end

end

TNew(1) = 500

TNew(nx) = 298

ENew(1)= K./(TNew(1)-TC(1))

ENew(nx) = K/(2*(TC(nx)-TNew(nx)))

T = TNew

time = time + dt

figure(1),clf

[hAx,hLine1,hLine2] = plotyy(x,T,x,ENew);

title(['Dielectric Constant and Temperature evolution after ' ,num2str(time/dt),' Seconds'])

xlabel('Distance (m)', 'FontSize',12)

ylabel(hAx(1),'Temperature(K)', 'FontSize',12) % left y-axis

ylim(hAx(1),[250,550])

ylabel(hAx(2),'Dielectric Constant(\epsilon)', 'FontSize', 12)

drawnow

frame=getframe(gcf)

writeVideo(v,frame)

end

close(v)

figure(2)

plot(T,ENew)

title('Dielectric Constant as a Function of Temperature of BaTiO3')

xlabel('Temperature (K)','FontSize', 12)

ylabel('Dielectric Constant (\epsilon)','FontSize', 12)