% Ssolutions to the heat flow equation 
% on one-dimensional rod length W

format compact;
W = 10;				% width of plate
Temp = 100;			% Constant temperature of rod, initially
tfin = 20;			% final time
alpha = 1;			% heat coef of the medium

NptsX=151;          % number of x pts
NptsT=151;          % number of t pts
Nf=200;              % number of Fourier terms
x=linspace(0,W,NptsX);
t=linspace(0,tfin,NptsT);
[X,T]=meshgrid(x,t);

fs=8;
figure(1)
clf

b=zeros(1,Nf);
U=zeros(NptsT,NptsX);

for n=1:Nf
    b(n)=(2*Temp/(n*pi))*(1-(-1)^n);  % Fourier coefficients
    Un=b(n)*exp(-(n*pi*alpha/W)^2*T).*sin(n*pi*X/W);  % Temperature(n)
    U=U+Un;
end
set(gca,'FontSize',[fs]);
surf(X,T,U);
shading interp
xlabel('x','Fontsize',fs); ylabel('t','Fontsize',fs); zlabel('u(x,t)','Fontsize',fs);axis tight
view([60 32])

figure(2)
clf

set(gca,'FontSize',[fs]);
surf(X,T,U);
shading interp
colormap(jet)
view([0 90])             %create 2D color map of temperature
xlabel('x','Fontsize',fs); ylabel('t','Fontsize',fs); zlabel('u(x,t)','Fontsize',fs);axis tight
colorbar
set(gca,'FontSize',[fs]);

plotfilm=1;
if(plotfilm==1)
 fs=16;
 figure(3)
 clf
 set(gca,'FontSize',[fs]);
 xl=min(min(x));
 xr=max(max(x));
 yt=ceil(max(max(U)));
 yb=min(min(U));
 ax=[xl xr yb yt];
 for j=1:length(U(:,1))
  plot(x,U(j,:),'LineWidth',3)
  axis(ax);
  xlabel('x','Fontsize',fs); ylabel('u(x,t)','Fontsize',fs); 
  mytext=['t=',num2str(t(j)),'/',num2str(tfin)];
  text(xl+0.1*(xr-xl),yt-0.1*(yt-yb),mytext,'Fontsize',fs)
  pause(0.005)
  if(j==1)
   fprintf('Press any key to continue')
   pause
  end
 end
end