% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% heatSln
% ____________________________________________________________________________
% 18.385 (Nonlinear Dynamical Systems).         MIT, Fall 1999. R. R. Rosales.
%
% Demo illustrating the solution of the heat equation T_t = T_xx with periodic
% boundary conditions on 0 <= x <= 2*pi and  initial conditions u(x, 0) = F(x)
% as given by the same options in the fourierSC script (see help fourierSC for
% more information). A spectral method is used to solve the equation.  
%
% At the end of the run:
%  x    (row array) contains the points in space.
%  fct  (""""""""") contains the initial values.
%  tt   (column array) contains the computed times.
%  uu   contains the solution, each row corresponding to an entry in tt.
%  dt   contains the separation between computed times.
%  t_f  contains the final time.  .... Thus: tt = (0:dt:t_f)'
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%========== Vector sizes and useful text variables: ==========================
M  = 512;
N  = 2*M;
ind = [0:M,(1-M):(-1)];
FN  = 'FontName';
HNB = 'HelveticaNarrowBold';
FS  = 'FontSize';
LW  = 'LineWidth';
MS  = 'MarkerSize';
Space = '                      ';
%
%========== Enter inputs needed and print on-screen info: ====================
%
fprintf('\n Solution to the heat equation T_t = T_xx.')
fprintf('\n The solution will be plotted for t=0:dt:t_f \n')
%
dt  = input([Space, 'Enter dt ......................... dt = ']);
z   = exp(-dt*ind.^2);
t_f = input([Space, 'Enter t_f ....................... t_f = ']);
%
%========== Construct array of x's ===========================================
%
dx = 2*pi/N;
x  = 0:dx:2*pi-dx;
%
%========== Construct function ===============================================
%
fprintf('\n NOW CHOOSE AN INITIAL CONDITION. \n')
FSoption;
%
%========= Plotting Ranges ===================================================
%
Max = max(fct);
Min = min(fct);
range = Max - Min;
Max = Max + 0.1*range;
Min = Min - 0.1*range;
%
%========= READY TO CALCULATE AND PLOT =======================================
%
fprintf('\n PRESS ANY KEY TO MOVE FORWARD IN TIME.\n')
u = fct; uu = u;
t=0;     tt = 0;
figure
whitebg('k')
plot(x, u, '-r', LW, 3)
title('Heat equation: initial condition.', FN, HNB, FS, 22)
xlabel('x', FN, HNB, FS, 18)
ylabel('Solution', FN, HNB, FS, 18)
grid on
axis([0 2*pi Min Max])
pause
%
while t <= t_f - 0.1*dt
   u = real(ifft(fft(u).*z)); uu = [uu; u];
   t = t+dt; tt = [tt; t];
   TEXT = num2str(t);
   plot(x, u, '-r', LW, 3)
   title(['Heat equation solution, time t = ', TEXT], FN, HNB, FS, 22)
   xlabel('x', FN, HNB, FS, 18)
   ylabel('Solution', FN, HNB, FS, 18)
   grid on
   axis([0 2*pi Min Max])
   pause
end
figure
whitebg('k')
mesh(x, tt, uu)
shading interp
colormap(hsv)
axis([0 2*pi 0 t_f Min Max])
grid on
title('Solution of the heat equation.', FN, HNB, FS, 22)
xlabel('Space', FN, HNB, FS, 18)
ylabel('Time', FN, HNB, FS, 18)
zlabel('Solution', FN, HNB, FS, 18)
fprintf('\n DONE. \n')
%  
clear MM        OPTION    dyt       smooth    z         t         MS
clear Space     FN        Max       TEXT      h         FS        Min
clear ans       h_fig     u         HNB       N         ind       LW
clear NN        dx        jn        yb        M         OPT       dyb
clear range     yt        BGC       CLR
%
% EOF