myfunction not working when applied on a meshgrid

I am absolute beginner in Matlab, and I have a function that solves the kepler equation via newton's method (iterations), that works 100% correctly when I manually input argument, but when I try to apply it on a meshgrid elements, it doesn't work, I tried multiple solutions but to no avail.

the function in kepler.m

function [E_n,i_iter] = kepler(M_an,ecc_1)

f=@(E) M_an - E + (ecc_1*sin(E));
df= @(E) -1 + (ecc_1*cos(E));

error=0;
E0=0;
i=0;
while error > (10^-6)
    E1= E0- (f(E0)/df(E0));
    error = abs(E1-E0);
    E0=E1;
    i=i+1;
end
E_n=E1;
i_iter=i;

the script file with the grid

x= 0:0.01:2*pi; %629 element.
y= 0:0.001:1; %1001 element. 
[X, Y]= meshgrid(x,y);
z=zeros(1001,629);  %% for X and Y are both displayed 1001x629 double in the workspace.
z=kepler(X,Y); % it is here where it doesnt work.

% tried 
% for i=0:1001
%   for j=0:629
%       z(i,j)=(kepler(X(i,j),Y(i,j)));
%   end
% end 

when I do as in the code above I have the message

Warning: Rank deficient, rank = 1, tol =5.574419e-12.

> In kepler (line 9)

In untitled2 (line 6)

when I try the commented code I get :

Index in position 1 is invalid. Array indices must be positive integers or logical values.

Error in untitled2 (line 11)

z(i,j)=(kepler(X(i,j),Y(i,j)));

Thank you in advance for your help.

Answer 1

[Torsten]

x= 0:0.01:2*pi; %629 element.
y= 0:0.001:1; %1001 element. 
[X, Y]= meshgrid(x,y);
z=zeros(1001,629);  %% for X and Y are both displayed 1001x629 double in the workspace.
for i=1:1001
    for j=1:629
        [z(i,j),~]=kepler(X(i,j),Y(i,j));
    end
end 

function [E_n,i_iter] = kepler(M_an,ecc_1)
  f = @(E) M_an - E + (ecc_1*sin(E));
  df = @(E) -1 + (ecc_1*cos(E));
  error = 1;
  E0 = 0;
  i = 0;
  while error > (10^-6)
    E1 = E0 - f(E0)/df(E0);
    error = abs(E1-E0);
    E0 = E1;
    i = i+1;
  end
  E_n = E1;
  i_iter = i;
end