function pde1

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

m = 0;
x = linspace(0,pi,200);   % spatial grid
t = linspace(0,3,151);     % output times

xa = linspace(0,pi,40);   % spatial grid
ntrunc=5;
a=[];
for n=1:ntrunc
%	a(n)=trapz(x,sin(n*x).*sin(x).^2)*2/pi;
%[n,a(n),(1-(-1)^n)*(2/n-1/(n+2)-1/(n-2))/2/pi]
	a(n)=trapz(x,sin((n-1/2)*x).*sin(x))*2/pi;
[n,a(n),(-1)^(n+1)*(1/(n+1/2)-1/(n-3/2))/pi]
end

sol = pdepe(m,@pdex1pde,@pdex1ic,@pdex1bc,x,t);
u = sol(:,:,1);

subplot('position',[0.1 0.42 0.88 0.52])
%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');
xlabel('Position x');
ylabel('Time t');

ns=[1 3 7 21 51 151];
subplot('position',[0.1 0.09 0.4 0.23])
plot(x,u(ns,:),'b');
xlabel('Position x');
ylabel('u(x,t)');
axis([0 pi -inf inf])
text(1.,2.7,'t=0, 0.04, 0.12, 0.4, 1, 3')

hold on
  for nt=ns
  ua=0*xa;
  for n=1:ntrunc
	  ua=ua+a(n)*sin((n-1/2)*xa)*exp(-(n-1/2)^2*t(nt));
  end
  plot(xa,ua,'r.')
  end
hold off

subplot('position',[0.58 0.09 0.4 0.23])
plot(t,u(:,end),t,u(:,100));
xlabel('Time, t');
ylabel('u(\pi,t) and u(\pi/2,t)');
axis([0 t(end) -inf inf])


hold on
  for tt=t(1:5:end)
  u0 = 0; u1=0;
  for n=1:ntrunc
	  u0=u0+a(n)*sin((n-1/2)*pi)*exp(-(n-1/2)^2*tt);
          u1=u1+a(n)*sin((n-1/2)*pi/2)*exp(-(n-1/2)^2*tt);
  end
  plot(tt,u0,'b.',tt,u1,'r.')
  end
hold off


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

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

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

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

function u0 = pdex1ic(x)
  u0 = sin(x);
    
% --------------------------------------------------------------------------

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

pl = ul;
ql = 0;
pr = 0;
qr = 1;
  
