function Heat1D
%
% SOLUTION OF THE 1D HEAT EQUATION $u_{tt} - D u_{xx} = 0$
% IN $(-L,L) \times(0,T)$ WITH ZERO DIRICHLET CONDITIONS
%

format short;
% format long;
close all;

global D IPLOT NX1 T U0 NT1 XNUMER HX NX

%
% THE DATA OF THE PROBLEM
%

T=10.;                          % FINAL TIME
Lmax=1.;                        % HALF-LENGTH

D=.1;                           % DIFFUSION COEFFICIENT
u0max=.1; sigma0=.1;        % DATA FOR INITIAL DATA

NX=50;                         % FOR SPATIAL APPROXIMATION
% NX=20;                         % FOR SPATIAL APPROXIMATION
NT=50;                         % NUMBER OF TIME STEPS
% NT=20;                         % NUMBER OF TIME STEPS

IPLOT=2;                       % FIGURE REPRESENTATION LEVEL
% IDISP=1;                       % PRINT LEVEL

%
% PREPARING COMPUTATIONS
%

HX=2.*Lmax/NX;                  % SPACE STEP
NX1=NX+1;
XNUMER=linspace(-Lmax,Lmax,NX1);% THE x_j

% HT=T/Ntt;

coef0=.5/sigma0;                
U0=u0max*exp(-coef0*XNUMER.*XNUMER); % INITIAL CONDITION

% PIS2=.5*pi;
% U0=u0max*(sin(PIS2*(XNUMER+1.)));

NT1=NT+1;

%

tnumer=linspace(0.,T,NT1);      % THE t_i
sol=pdepe(0,@PDEfun,@ICfun,@BCfun,XNUMER,tnumer); % PDEPE

% ATTENTION: SOL(j,k) IS AN APPROXIMATION TO THE SOLUTION AT
%            (t,x)=(tnumer(j),XNUMER(k))
    
if (IPLOT > 0)
    
    figure(1);
    view(-70.,20.);             % ANGLES: AZIMUTH, ELEVATION
    surf(XNUMER,tnumer,sol);
    xlabel('x: SPACE');
    ylabel('t: TIME');
    zlabel('u: TEMPERATURE');
    
    figure(2);
    maxx=max(sol');
    plot(tnumer,maxx,'LineWidth',1.2);
    legend('SPATIAL MAXIMUM');
    axis tight;
    xlabel('TIMES');
    ylabel('MAX SOL');
    title(['INITIAL AND FINAL MAXIMAL VALUES: ',num2str(maxx(1)),',  ',...
        num2str(maxx(end))]);
    shg;
    
end

% format long;

end

function [c,f,s]=PDEfun(~,~,~,Du) % INDICATING C, F, S

% [C,F,S]=PDEFUN(X,T,U,DUDX)

global D IPLOT NX1 T U0 NT1 XNUMER HX NX

c=1.;
f=D*Du;
s=0.;
end

function u=ICfun(x) % INDICATING U AT TIME 0

%

global D IPLOT NX1 T U0 NT1 XNUMER HX NX

u=interp1(XNUMER,U0,x);
end

function [pl,ql,pr,qr]=BCfun(~,ul,~,ur,~) % INDICATING BOUNDARY CONDITIONS

% [PL,QL,PR,QR]=BCFUN(XL,UL,XR,UR,T)

pl=ul; ql=0.;
pr=ur; qr=0.;
end
