function plotu

%% this plots the solution u(x,t) for the hot-spot on an infinite rod (section 1,
%% fourtran.pdf)

%% to use this for other u(x,t), change the function uxt below
%% you may also want to change L, the length of the interval

close all; % close all open figures

% setup
L=5;
x=linspace(-L,L,100);

%% plot initial condition u(x,0) vs. x
figure(1);
hdl0=plot([-L,-1/2],[0,0],'k',[-1/2 -1/2],[0,1],'k--',[-1/2,1/2],[1,1],'k',[1/2 1/2],[0,1],'k--',[1/2,L],[0,0],'k',...
    -1/2,0,'ko',-1/2,1,'ko',1/2,0,'ko',1/2,1,'ko'); hold on;

%% plot u(x,t0) for t0=1/4,1
hdl1=plot(x,uxt(x,1/4),'b');
hdl2=plot(x,uxt(x,1),'r');

% formatting the figure
set(hdl0,'linewidth',1);
set(hdl1,'linewidth',1);
set(hdl2,'linewidth',1);
set(gca,'fontname','times','fontsize',14,'xlim',[-L,L],'ylim',[0,1.1]);
xlabel('x','fontname','times','fontsize',14);
ylabel('u(x,t_0)/u_0 (\epsilon=1)','fontname','times','fontsize',14);

% plot u(x0,t) vs. time t for different fixed x0
Tf=5; % max time to plot
t=linspace(0.001,Tf,100);
figure(2),plot(t,uxt(0,t),'k',t,uxt(1/4,t),'k--',t,uxt(3/4,t),'k-.');
% formatting the figure
set(gca,'fontname','times','fontsize',14);
xlabel('\kappa t','fontname','times','fontsize',14);
ylabel('u(x_0,t)/u_0 (\epsilon=1)','fontname','times','fontsize',14);

% contour plot of u(x,t)
cvals=[0.75,0.5,0.25,0.2,0.1,0.05]; % values of the contours to plot, i.e. u(x,t)=const=cvals
[xx,tt]=meshgrid(x,t);
uxt3D=uxt(xx,tt);
figure(3);
[cs,hdl]=contour(xx,tt,uxt3D,cvals); 
% numbering the contours
clabel(cs,hdl,'fontsize',15,'rotation',0,'labelspacing',1024)
% formatting the figure
set(gca,'xlim',[-L,L],'ylim',[0,Tf],'fontname','times','fontsize',14);
xlabel('x','fontname','times','fontsize',14); 
ylabel('\kappa t','fontname','times','fontsize',14); 


% 3D plot of u(x,t)
figure(4),mesh(xx,tt,uxt3D);
set(gca,'xlim',[-L,L],'ylim',[0,Tf],'fontname','times','fontsize',14,'view',[150,40]);
xlabel('x','fontname','times','fontsize',14); 
ylabel('\kappa t','fontname','times','fontsize',14); 
zlabel('u(x,t)','fontname','times','fontsize',14); 

%------------------------------------------
% your solution u(x,t) - change this for different problems
function u=uxt(x,t)

% epsilon=1, kappa=1
u=1/2*(erf( (x+1/2)./(2*sqrt(t)) ) - erf( (x-1/2)./(2*sqrt(t))));