function pd23a

' Solves the PDE using an in-built function PDEPE'

m = 0;
x = linspace(0,pi/2,200);   % spatial grid
t = linspace(0,1,201).^2*4; % time grid
ts = [0 .01 0.09 0.36 1 2]; % special output times
t = unique(sort([ts t]));   % all output times
ns = [];
for tt = ts, ns = [ns find(t==tt)]; end

sol = pdepe(m,@pdex1pde,@pdex1ic,@pdex1bc,x,t);    % solve PDE
u = sol(:,:,1);            % the solution...

% Everything else is output

subplot('position',[0.1 0.48 0.88 0.46])
%imagesc(x,t,u),axis xy
%contour(x,t,u,20)
contourf(x,t,u,20)
  colorbar
%pcolor(x,t,u),shading interp

title('Numerical solution computed with 200 mesh points','fontsize',12);
xlabel('Position x','fontsize',12);
ylabel('Time t','fontsize',12);

set(gca,'yscal','log')

subplot('position',[0.1 0.09 0.4 0.29])
plot(x,u(ns,:),'b');
xlabel('Position x','fontsize',12);
ylabel('u(x,t)','fontsize',12);
axis([0 pi/2 0 1.02])
text(.2,1.1,'t=0, 0.01, 0.09, 0.36, 1, 2','fontsize',12)


subplot('position',[0.58 0.09 0.4 0.29])
plot(t.^.5,4/pi*exp(-t),'r--','linewidth',1);
hold on
plot(t.^.5,u(:,1),'b')
plot(t.^.5,u(:,133)*2,'color',[0 .5 0])
hold off
xlabel('(Time,t)^{1/2}','fontsize',12);
ylabel('u(0,t) and 2u(\pi/3,t)','fontsize',12);
axis([0 t(end).^.5 0 inf])
text(.03,.8,'\uparrow u(0,t)','fontsize',13,'color','b')
text(.3,1.71,'\leftarrow 2u(\pi/3,t)','fontsize',13,'color',[ 0 .5 0])


% --------------------------------------------------------------------------

function [c,f,s] = pdex1pde(x,t,u,DuDx)

	% The PDE

c = 1;
f = DuDx;
s = 0;
end

% --------------------------------------------------------------------------

function u0 = pdex1ic(x) 

	% The initial condition

  u0 = 1;
end

% --------------------------------------------------------------------------

function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t)

	% The boundary conditions
  
pl = 0;
ql = 1;
pr = ur;
qr = 0;
end

end
