x = linspace(-2, 3, 501);
y = x.^2 - 1;
ind = find(y>2);
plot(x, y, 'b-', x(ind), y(ind), 'rx');