function q=DblSimpson(f,a,A,b,B,m,n)

if(m==1 && n==1)           %辛普森公式

q=((B-b)*(A-a)/9)*(subs(sym(f),findsym(sym(f)),{a,b})+...

subs(sym(f),findsym(sym(f)),{a,B})+...

subs(sym(f),findsym(sym(f)),{A,b})+...

subs(sym(f),findsym(sym(f)),{A,B})+...

4*subs(sym(f),findsym(sym(f)),{(A-a)/2,b})+...

4*subs(sym(f),findsym(sym(f)),{(A-a)/2,B})+...

4*subs(sym(f),findsym(sym(f)),{a,(B-b)/2})+...

4*subs(sym(f),findsym(sym(f)),{A,(B-b)/2})+...

16*subs(sym(f),findsym(sym(f)),{(A-a)/2,(B-b)/2}));

else                       %复合辛普森公式

q=0;

for i=0:n-1

for j=0:m-1

x=a+2*i*(A-a)/2/n;

y=b+2*j*(B-b)/2/m;

x1=a+(2*i+1)*(A-a)/2/n;

y1=b+(2*j+1)*(B-b)/2/m;

x2=a+2*(i+1)*(A-a)/2/n;

y2=b+2*(j+1)*(B-b)/2/m;

q=q+subs(sym(f),findsym(sym(f)),{x,y})+...

subs(sym(f),findsym(sym(f)),{x,y2})+...

subs(sym(f),findsym(sym(f)),{x2,y})+...

subs(sym(f),findsym(sym(f)),{x2,y2})+...

4*subs(sym(f),findsym(sym(f)),{x,y1})+...

4*subs(sym(f),findsym(sym(f)),{x2,y1})+...

4*subs(sym(f),findsym(sym(f)),{x1,y})+...

4*subs(sym(f),findsym(sym(f)),{x1,y2})+...

16*subs(sym(f),findsym(sym(f)),{x1,y1});

end

end

end

q=((B-b)*(A-a)/36/m/n)*q;