
function compute_integrals

f=@(x,y) exp(-x.^2-y.^2);
% f=@(x,y) (x.^2+y.^2).^(1/2);
% f=@(x,y) ((x-0.5).^2+(y-0.5).^2).^(1/2);
F=chebfun2(f,[0 1 0 1]);
I=sum2(F);
Q = integral2(f,0,1,0,1,'AbsTol',10^(-15));

format long e
I
Q
I-Q

plot(F)