Tricky. Suppose that your fit is
y = p(1)*x^3 + p(2)*x^2 + p(3)*x+p(4);
You'll need to ensure that the first derivative
yp = @(x) 3*p(1)*x^2 + 2*p(2)*x^2 + p(3);
is positive across the entire interval [x1,x2]. Certainly
yp(x1) >= 0 && yp(x2) >= 0
must be satisfied. If the slope doesn't have a minimum in the interval, you're done. Taking the derivative of yp shows that an extremum occurs at
(for this analysis it doesn't really matter if it is a minimum or maximum). If
the extremum is outside the interval. Otherwise, you need to calculate the minimum and ensure it is nonnegative:
Putting the constraints together,
-yp(x1) <= 0
-yp(x2) <= 0
(x1 <= x0 && x0 <= x1)*(p(2)^2-3*p(1)*p(3)) <= 0
You could use this to create a constraint function c(p) that you can use in fmincon (you would be minimizing the sum of squares).
On the other hand, if your data are monotonically increasing, you could discard the cubic fit and interpolate using pchip.