How can I get Octave to correctly compute a double integral of a piecewise function?

Computational Science People:

The title is the question. I am trying to numerically compute a certain integral over a square in Octave (an open-source Matlab clone). I'm getting the wrong answer. To show that something really is going wrong, I've created the following toy problem: Let $f = \chi_{(-\infty, 1)}$ ($f(t) = 1$ if $t < 1$, otherwise $0$). Then $\int_0^3 \int_0^3 f(t) f(t')\,dt'\,dt = 1$. I wrote an Octave script to compute this:

    junk=10;  %this is here just so Octave doesn't complain about the file name
function out =R1(t)
if t < 1
out = 1;
else
out =0;
end%if
end%function

dIF = @(t1,t2) ones(size(t1))* R1(t1(1))*R1(t2);

dblquad(dIF, 0, 3, 0, 3)  %should be 1.  Octave incorrectly gives 3.

tx=0:0.1:3;
ty=0:0.1:3;
[xx, yy] = meshgrid (tx, ty);
z=zeros(31,31);

for i=1:31
for j=1:31
z(i,j)=dIF(xx(i,j),yy(i,j));
end%for j
end%for i
mesh(xx,yy,z) %graph is correct

%end script


Octave returns $3$ for the value of the double integral, and it should be $1$. As you can see, the script plots the integrand, which is, as intended, $1$ if $x$ and $y$ are both between $0$ and $1$, and $0$ otherwise. I defined the function dIF the way I did because in my actual problem, I need to pass additional parameters to the function I'm integrating, and this is the only way I know how to do it.

In my actual problem, I'll probably just cut my double integral into nine pieces, not requiring piecewise functions, and add them. But I'd like to know if there's a better way to do it, in case I face a similar problem in the future.

I have not attempted this in Matlab, only Octave. I tested it just now with a copy of Octave I downloaded today from what I am quite sure is a reputable source.

Your R1 function should be vector valued (it appears as though it is scalar valued). For example, this appears to work:

function out = R1(t)
out = t < 1;
end

dIF = @(x,y) R1(x) .* R1(y);