Speeding up pdepe calculation

I have a pdepe bit of code I want to speed up if I can. I have a (typically) 400x800 array as a source term, the code I'm currently using is:

function sol=T_sol(t,r,X,T_bulk,T_bdy)

m=1; %Sets the geometry to cylindrical
global theta kappa h Q W;
theta=X(1); kappa=X(2); h=X(3); %Parameters used
[t_grid,r_grid] = meshgrid(t,r);
Q = @(z) interp1(t,T_bdy,z);
W=@(tq,rq) interp2(t_grid,r_grid,T_bulk',tq,rq)';
sol=pdepe(m,@pdefun,@icfun,@bcfun,r,t); %Solve the PDE
end

%Defining PDE
function [c,f,s] = pdefun(r,t,u,DuDx)
global theta kappa W;
c = theta;
s = W(t,r);
f = kappa*DuDx;
end

%Initial conditions
function u0 = icfun(r)
u0 = 0;
end

%Boundary conditions
function [pl,ql,pr,qr] = bcfun(xl,ul,xr,ur,t)
global h Q;
pl = 0;
ql = 0;
pr = h*ur+Q(t);
qr = 1;
end

Currently it's taking about 5 minutes to run. Is there anyway I can speed it up?

[Matthew Hunt]

There was one run which took 2 and a half hours!!!