clear all
%implicit 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=10; %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=300; %number of time steps to do
time(1)=0;
tmid(1)=0; %initial centerline temperature
temperature=zeros(ndiv-1); %fill temperature vector with zeros
a=(1+2*alpha)*eye(ndiv-1)-alpha*diag(ones(ndiv-2,1)',1)-alpha*diag(ones(ndiv-2,1)',-1);
for i=1:ntimes
    b=temperature;
    b(1)=b(1)+alpha*tboundary;
    b(ndiv-1)=b(ndiv-1)+alpha*tboundary;
    time(i+1)=i*dt;
    temperature=a\b;
    tmid(i+1)=temperature(ndiv/2);
end
plot(time,tmid)
