Notice that j is empty, then try this:
X = -1.41:0.01:1.41;
[~,idx] = min(abs(X-1));
disp(X(idx)-1)
This being floating point, you're not getting exact increments of 
. You could scale everything so x, y are 
and then find 
and 
. Or you could do it like this:
. You could scale everything so x, y are and then find and . Or you could do it like this:z = @(x,y) sqrt(4-x.^2-y.^2);
X = -1.41:0.01:1.41;
[x,y] = meshgrid(X);
[fx,fy] = gradient(z(x,y),0.01);
[~,j] = min(abs(X-1));
fx0 = fx(j,j);
fy0 = fy(j,j);
z1 = @(x,y) z(1,1) + fx0*(x-1) + fy0*(y-1);
As a side note, you don't need 
. You can use j directly as logical indexing.
. You can use j directly as logical indexing.
