It's largely correct, except that your function distsq is not differentiable at x=0. So, if there's a chance the solution might lie there (but I think it's impossible if P(n) and at least one other P(i) are greater than zero), then I would make a transformation to get rid of the non-differentiability. In this case, this could be,
distsq = @(x) sum((x-P(1:end-1)).^2) + (f(x).^2-P(end).^2).^2;
Note however that for the specific f in your example, the transformation turns the problem into a linear least squares problem, so that you can use lsqlin instead of fmincon,
C=[speye(n-1);ones(1,n-1)];
d=P; d(end)=d(end)^2;
[x,fval] = lsqlin(C,d,[],[],[],[],[0 0],[1 1]);
This also has the advantage that lsqlin is globally convergent and doesn't require an initial guess.
