You need to build what is called a delayed differential equation. I was about to explain how to do, but then I found this wonderful tutorial to do just that. Example 1 is basically what you need.
The only extra caveat is that you should incorporate dx/dt and dy/dt into the same set of differential equation
Let me know if you need more help
Edit: keep it in one file
function dYdt = ddefun(t,Y,Z)
% assume Y = [x;y]
x = Y(1:n); % 2n is the size of Y. this step is unnecessary ...
y = Y(n+1:2*n); % but helps visualize what is happening
ytau = Z(:,1);
dYdt(1:n) = -c*x.*y;
dYdt(n+1:2*n) = a*dot(x,y) + x.*ytau + b
end