clear all
%explicit parabolic PDE
%
k=100; %W/m^2
rho=1000; %kg/m^3
cp=500; %J/kg-K
diffus=k/rho/cp; %thermal diffusivity
length=1; %m
dt=0.2; %s
tboundary=100; %boundary temperature (C)
ndiv=100; %number of spatial divisions
h=length/ndiv; %mesh spacing
%
alpha=diffus*dt/h^2 %have this print to check stability (<0.5)
%
ntimes=15000; %number of time steps to do
time(1)=0;
tmid(1)=0; %initial centerline temperature
told=zeros(ndiv+1); %fill temperature vector with zeros
for i=1:ntimes
    time(i+1)=i*dt;
    temperature(1)=tboundary;
    temperature(ndiv+1)=tboundary;
    for j=2:ndiv
        temperature(j)=(1-2*alpha)*told(j)+alpha*(told(j-1)+told(j+1));
    end
    told=temperature;
    tmid(i+1)=temperature(ndiv/2+1);
end
plot(time,tmid)
