If you expand out your expression, you'll see that you have a polynomial up to x^3, so there are four coefficients including the intercept. But you have specified five coefficients. Try this instead:
modelFun = @(p,x) p(1)*(x-p(4)) + p(2)*(x-p(4)).^2 + p(3)*(x-p(4)).^3;
p0 = [1 1 1 500];
Then use p0 as your initial values. You may find that both nlinfit and lsqcurvefit will have an easier time of it if you avoid trying to include the fifth (overdetermined) parameter. Then, if you like that, you can notice that you get the same results this way:
poly = polyfit(x,y,3);
line(x,polyval(poly,x),'color','c')
Your function is equivalent to a third-order polynomial, and it can be fit by linear least squares with no nonlinear fitting required.
