%               Function simp2
%
%  Two dimensional Simpson's rule over a rectangle.
%  The call is simp2(f,corners) when f is given as an inline function,
%  and simp2('f',corners) when f is given in an mfile. 
%      corners = [a b c d] is a vector of corner coordinates of the
%  rectangle with corners (a,c), (b,c), (b,d), and (a,d). We assume
%  a< b and c < d.
%     After the call the user is asked to enter the number n of 
% subdivisions in the x direction, and the number m of the subdvision
% in the y direction.  n and m must be even. 

function out = simp2(f, corners)

 a = corners(1); b = corners(2); c = corners(3); d = corners(4);
 disp('  ')
 disp('the number of subivisions n and m in each direction must be even ')
 subdiv = input('enter the number of subdivisions in x and y direction as [n m] ')
 n = subdiv(1); m = subdiv(2);

 x = linspace(a,b, n+1); y = linspace(c,d,m+1);
 [X,Y] = meshgrid(x,y);
 
 svecx = 2*ones(size(x));
 svecx(2:2:n) = 4*ones(1,n/2);
 svecx(1) = 1; svecx(n+1) = 1;

 svecy = 2*ones(size(y));
 svecy(2:2:m) = 4*ones(1,m/2);
 svecy(1) = 1; svecy(m+1) = 1;

 S = svecy'*svecx;
 T = S.*feval(f,X,Y);
 disp('Approximate value of the integral using Simpsons rule ')
 out = sum(sum(T))*(b-a)*(d-c)/(9*m*n);