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

format compact;
W = 10;				% width of plate
Temp = 100;			% temperature of RHS of slab
alpha = 1;			% heat coef of the medium

NptsT=101;
Nplots=500;Nrows=5;plotsols=0;tfin=40;NptsX=401;
%Nplots=6;Nrows=2;plotsols=1;tfin=5;NptsX=201;
Ncols=Nplots/Nrows;
Nskip=1;			% use every Nskip coefs
x=linspace(0,W,NptsX);
t=linspace(0,tfin,NptsT);
[X,T]=meshgrid(x,t);


fs=8;

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);
%  b(jj) = -200*((-1)^(jj)+1)/(jj*pi);
end

b(1:Nskip:Nplots)

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

if(plotsols==1)
  figure(3);clf
end

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)
   figure(3);
   set(gca,'FontSize',[fs]);
   subplot(2*Nrows,Ncols,jj+Ncols*floor((jj-1)/Ncols))
   set(gca,'FontSize',[fs]);
   surf(X,T,Un);
   shading interp
   %colormap(gray)
   zl = ['z(b_{',num2str(jj),'})'];
   xlabel('x','Fontsize',fs); ylabel('t','Fontsize',fs); zlabel(zl,'Fontsize',fs);axis tight
   view([12 46])
   view([90 0])
   if(jj==1)
     ax=axis;
     ax(5)=-ax(6)/3; 
   end
   axis(ax);
   axis tight;
   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)
   drawnow
  end
end

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

% Now plot the Fourier combo:
fs=16;
figure(4);clf
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),'%3.1f'),'/',num2str(tfin)];
  text(xl+0.1*(xr-xl),yt-0.1*(yt-yb),mytext,'Fontsize',fs)
  %pause(0.001)
  drawnow
  if(jj==1)
   fprintf('Press any key to continue\n')
   pause
   figure(5)
  end
 end
end



