% Generating solutions to the heat flow equation 
% in an infinite slab of thickness W

format compact;
W = 10;				% width of plate
Temp = 200;			% temperature of RHS of slab
tfin = 20;			% final time
alpha = 1;			% heat coef of the medium

NptsX=151;
NptsT=151;
Nplots=500; Nrows=5;			% number of columns
Nplots=540; Nrows=6;			% number of columns
plotsols=0;
Ncols=Nplots/Nrows;
Nskip=2;			% use every Nskip coefs
x=linspace(0,W,NptsX);
t=linspace(0,tfin,NptsT);
[X,T]=meshgrid(x,t);


fs=8;
figure(3)
clf

U0t = Temp/W*x;			% initial theoretical steady state distribution
b=zeros(1,Nplots);

% Fourier coefs for U0(x) temperature
for jj=1:Nskip:Nplots
  b(jj) = (2*Temp)/(jj*pi)*(-1)^(jj-1); 		% Triangular IC
%  b(jj) = -200*((-1)^(jj)+1)/(jj*pi);
  b(jj) = Temp*(1-cos(jj*pi))/(jj*pi);			% Square IC
end

b(1:Nskip:Nplots)

U=zeros(NptsT,NptsX);			% initial steady state distribution (Fourier)

for jj=1:Nskip:Nplots		% build Fourier series
  Un=b(jj)*exp(-(jj*pi*alpha/W)^2*T).*sin(jj*pi*X/W);
  U=U+Un;
  if(plotsols==1)
   subplot(2*Nrows,Ncols,jj+Ncols*floor((jj-1)/Ncols))
   set(gca,'FontSize',[fs]);
   surf(X,T,Un);
   shading interp
   %colormap(gray)
   zl = ['B_{',num2str(jj),'} \phi_{',num2str(jj),'}(x)'];
   xlabel('x','Fontsize',fs); ylabel('t','Fontsize',fs); zlabel(zl,'Fontsize',fs);axis tight
   view([12 46])
   if(jj==1)
     ax=axis;
     ax(5)=-ax(6)/3; 
   end
   axis(ax);
   set(gca,'FontSize',[fs]);
  
   subplot(2*Nrows,Ncols,jj+Ncols*floor((jj-1)/Ncols)+Ncols)
   set(gca,'FontSize',[fs]);
   surf(X,T,Un);
   shading interp
   colormap(jet)
   view([0 90])
   xlabel('x','Fontsize',fs); ylabel('t','Fontsize',fs); zlabel('z','Fontsize',fs);axis tight
   set(gca,'FontSize',[fs]);
   %fprintf('Press any key to continue...\n');
   %pause
   pause(0.0001)
  end
end

%Uinf=100-100*X/W;
%U=U+Uinf;

% Now plot the Fourier combo:
fs=16;
figure(4)
clf

  subplot(2,1,1)
  set(gca,'FontSize',[fs]);
  surf(X,T,U);
  shading interp
  %colormap(gray)
  xlabel('x','Fontsize',fs); ylabel('t','Fontsize',fs); zlabel('z','Fontsize',fs);axis tight
  view([60 32])
  set(gca,'FontSize',[fs]);

  subplot(2,1,2)
  set(gca,'FontSize',[fs]);
  surf(X,T,U);
  shading interp
  colormap(jet)
  view([0 90])
  xlabel('x','Fontsize',fs); ylabel('t','Fontsize',fs); zlabel('z','Fontsize',fs);axis tight
  colorbar
  set(gca,'FontSize',[fs]);

plotfilm=1;
if(plotfilm==1)
 fs=16;
 figure(5)
 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 jj=1:length(U(:,1))
  plot(x,U(jj,:),'LineWidth',3)
  axis(ax);
  xlabel('x','Fontsize',fs); ylabel('u(x,t)','Fontsize',fs); 
  mytext=['t=',num2str(t(jj)),'/',num2str(tfin)];
  text(xl+0.1*(xr-xl),yt-0.1*(yt-yb),mytext,'Fontsize',fs)
  pause(0.001)
  if(jj==1)
   fprintf('Press any key to continue')
   pause
  end
 end
end