define_variable(dblint_y,10,fixnum,"number of y divisions");
define_variable(dblint_x,10,fixnum,"number of x divisions");
dblint(f,c,d,a,b):=
     (mode_declare([function(f,c,d),a,b],float),
		  block([m2,n2,h,j1,j2,j3,x,dox,cox,hx,k1,k2,k3,y,z,l],
			mode_declare([m2,n2,h,j1,j2,j3],float),
				n2:.5/dblint_x,
	m2:.5/dblint_y,
				h:(b-a)*n2,
				j1:0.,j2:0.,j3:0.,
				for i:0 thru 2*dblint_x step 1 do
				(mode_declare(i,fixnum,[x,dox,cox,hx,k1,k2,k3],float),
					     x:a+float(i)*h,
					     dox:apply(d,[x]),cox:apply(c,[x]),
					     hx:(dox-cox)*m2,
					     k1:apply(f,[x,cox])+
                                                apply(f,[x,dox]),
					     k2:0.,k3:0.,
				    for j:1 thru 2*dblint_y-1 step 1 do
				    (mode_declare(j,fixnum,[y,z,l],float),
					     y:cox+float(j)*hx,
					     z:apply(f,[x,y]),
					     if evenp(j) then k2:k2+z else k3:k3+z),
				    l:(k1+2.*k2+4.*k3)*hx/3.,
				    if (i=0 or i=2*dblint_x) then j1:j1+l else
				       if evenp(i) then j2:j2+l else j3:j3+l),
				return((j1+2.*j2+4.*j3)*h/3.)));