function z = simp( x, y)

if nargin<2
	y=x;
end;
[n m ] = size( y);
if n == 1
	y = y(:);
	n = m;
end;
if nargin < 2
	x = 1:n;
end;
x = x(:);
if length( x) ~= n
	error( 'Inputs must be of same length!');
end;
if rem( n, 2) == 0
	error( 'Uneven numer of intervals!');
end;

c = [ 1 4 1];
dx = diff( x);
z = 0;
for i = 1: 2: (n-2)
	if abs( dx(i) - dx(i+1)) > 0.00001
		error( 'Two sequential intervals not equal!');
	end
	z = z + c* [ y( i, :); y( i+1, :); y( i+2, :)] * dx(i);
end;
z = z/3;