MATLAB integral3 Error

Here is my Code, dexact_v3 is not working with some issue. Your help will highly appreciated.

t=0;L = 1; hx=0.01;
a=0;b=L;
x=a:hx:b;y=x;
My_result_W=wexact_Patrick(t,x,y);
%My_result_d=dexact_v3(t,x,y);




function [d1]=dexact_v3(t,x,y);
t=0;
L = 1; % Length of the box
hx=0.01;
a=0;
b=L;
x=a:hx:b;
y=x;
Ntrunc = 10; % Number of summations (high frequency oscillations can be neglected)
T = 10; % Duration of the time interval
p=1;
q=1;
Ui = @(x,y) (cos(pi.*x).*cos(pi.*y)); % An initial distribution for u, i.e. the u(x,y,0)
Vi = @(x,y) (sin(pi.*x).*sin(pi.*y)); % An initial distribution for v, i.e. the v(x,y,0)
di = @(x,y) Ui(x,y) - Vi(x,y); % Initial condition for w, i.e. w(x,y,0)

% Computation of the b'-coefficients (denoted by bd(p,q))

for p=1:Ntrunc
	for q =1:Ntrunc
		u_integrand = @(x,y) 4*di(x,y).*sin(p.*x).*cos(q.*y)/L.^2; 
		m_bd(p,q)= integral2(@(x,y)u_integrand(x,y),0,L,0,L); 
    end
end 
% Zeroth order of d(x,y,t) without the convolution integral term, this is denoted by dzeroh
dzeroh = 1; % Constant term
for p=1: Ntrunc
for q= 1:Ntrunc
			r = p.^2 + q.^2 + 1;
		dzeroh = @(x,y,t) (dzeroh(x,y,t) + m_bd(p,q).*exp(-r.*t).*sin(p.*x).*cos(q.*y)); 
		q = q+1;
	end
	p = p+1;
end 
hx = 0.1; % discretization of the x-interval
hy = 0.1; % " " " y- "
ht = 0.1; % " " " z- "
xc = a:hx:L;
yc = xc;
tc = 0:ht:T;
G= @(x,y,t,xc,yc,tc) ( sqrt(pi/(t-tc)).*e^(-(t-tc)).*exp(-((x-xc).^2+(y-yc).^2)/(4*(t-tc))) ); % The Green function
integrand = @(x,y,t,xc,yc,tc) (-1.*((wexact_Patrick(tc,xc,tc)+1).*(wexact_Patrick(tc,xc,yc)-1).^2).*G(x,y,t,xc,yc,tc)/4); 
dzeroadd = integral3(@(xc,yc,tc)integrand(xc,yc,tc),0,L,0,L,0,T); % Convolution term
% % Add the convolution integral term to the homogeneous solution (the solution without the nonlinear term)
dzero = dzeroh + dzeroadd; % It holds: dzero = 1 + d'_0

% Here the complete solution for d is computed

integrand = @(x,y,t,xc,yc,tc) (-1.*(wexact_Patrick(xc,yc,tc)+dzero(xc,yc,tc)).*(wexact_Patrick(xc,yc,tc)-dzero(xc,yc,tc)).^2.*G(x,y,t,xc,yc,tc)/4); 
dadd = integral3(@(xc,yc,tc) integrand,0,L,0,L,0,T); % Another convolution term to be added

% The final result for d1 is 
d1 = dzeroh + dadd;  
 return
%%%%%%%%%%%%%%%%%%%%%%% Working Code which I used in above Code%%%%%%%%%
function [W1]=wexact_Patrick(t,x,y);
L=1;
Ntrunc = 10;
p=1;
q=1;
Ui = @(x,y)  sin(pi*x).*sin(pi*y); % An initial distribution for u, i.e. the u(x,y,0)
Vi = @(x,y)  sin(pi*x).*sin(pi*y); % An initial distribution for v, i.e. the v(x,y,0)
Wi = @(x,y) Ui(x,y) + Vi(x,y); % Initial condition for w, i.e. w(x,y,0)
% Computation of the b-coefficients
for  p = 1:Ntrunc
for q = 1:Ntrunc
        my_integrand = @(x,y)  4*Wi(x,y).*sin(p.*x).*cos(q.*y)/L^2; 
        bd(p,q) = integral2(@(x,y)my_integrand(x,y),0,L,0,L);
    end
end
% The analytical solution for w(x,y,t) and first order of d(x,y,t)
W1 = 1; %Constant term
for p =1: Ntrunc
for q =1: Ntrunc
		r = p^2 + q^2 + 1;
		W1 = W1 + bd(p,q).*exp(-r.*t).*sin(p.*x).*cos(q.*y); 
    end
end 
return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%