How can i handle division by zero inside the symsum function?

I'm trying to model electric potential energy of 3 proton arranged in A(-1,0), B(0,0), C(0,1), graphically is(the next plot was made manually in geogebra 3d):

Asumming k=1, q1=1, q2=2 , mathematically can been described as:

The thing is that i want to generalize this into a summation function in matlab for n charges separated in the x-axis by a distance of 1 unit, mathematically it will be:

So i try the next code in matlab going form a=-1 to b=1:

[x,y]=meshgrid(-10:1:10);

syms k;
z=symsum((1./((x-k).^2+y.^2)),k,-1,1);

%Converts sym data type from the variable z to double data type for 3d plotting
doubleN = double(z);
surf(doubleN)

Next, it appears me the next code:

Answer 1

[Matt J]

Why do things symbolically at all?

[x,y,k]=meshgrid(-10:1:10,-10:1:10, -1:+1);

z= sum(  1./((x-k).^2+y.^2) , 3);

surf(x(:,:,1),y(:,:,1),z);

Answer 2

[Walter Roberson]

z=symsum((1./((x-k).^2+y.^2)),k,-1,1);

That is not worth doing with symsum. It is only 3 terms; just write them out, or use a vectorized calculation.

For example

dim = max(ndims(x),ndims(y)) + 1;
k = reshape(-1:1:1, 1, 1, dim);
z = sum(1./((x-k).^2+y.^2), dim);

In the case that y == 0 then Yes you get a division by 0 at each place one of the x exactly equal one of the k values -- but those will not stop the calculation (you will get +/- inf or nan depending on the exact expression.)