Use logical indexing.
function toint = toint(p, r, R, dRdt, mu)
c = 29979245800; %cm/s
gamma = 1/sqrt(2*(1 - (dRdt/c).^2));
limitmin = R - deltR(R, dRdt);
limitmax = R;
minsize = size(limitmin)
maxsize = size(limitmax)
%Preallocate the array toint as zeros, so that you only have to assign
%values to indices corresponding to the condition
toint = zeros(size(r));
%Indices corresponding to the condition
idx = (r > limitmin) & (r < limitmax);
%Assuming r is the only non-scalar array here, and rest of the variables are scalars
%If any other variables are non-scalar, use the index idx for them as well
toint = r(idx).^2./((gamma.^(p-1)).*(1-getbetastar(dRdt/c).*mu).^((p+3)/2));
end
