function [u,x,h] = SolvePDE(n, t, a, d, method)

   global exponent;
   
   h=1/n; % Grid spacing
   x=([1:n]-0.5)*h; % Cell centers
   IC = @(x) sin(pi*x).^exponent;
   u0 = FVQuad(IC,x,h);
   
   % For some silly reason MATLAB's odesolve switches the output format depending
   % on whether the vector of times has two or more elements (!?!)
   % So we add one more point here
   [T,W] = SolveODE(n, u0, [0 t/2 t], h, a, d, method);
   u=W(3,:);

end
