Plot surface in region with diffrent curvature on each boundary (square and circle), with acceptable resolution

I want to create a surface plot in a region that its outside boundary is a square and its inside a circle (black area of the following image):

This is my code:

x = -5:0.1:5;
y = -5:0.1:5;
A = 1;
d = 1;

[X, Y] = meshgrid(x,y);

R = sqrt(X.^2 + Y.^2);
Z = A ./ R;
REGION = X.^2 + Y.^2 >= d^2;
Z(~REGION) = NaN;

surface(X,Y,Z);
shading interp

If I choose cartesian coordinates system for the plot and then remove the inner circle with

Z(~REGION) = NaN;

Then I get an effect that looks like low resolution on the inner boundary:

If I choose the polar coordinates system (use instead of , so that the inner boundary appears alright) I have to either remove everything that is outside from , with the same command as above, which still produces the low resolution effect but on the square boundary, or plot beyond , and then focus the plot camera to the desired region (this option produces the desired effect). Is there an elegant way to plot surface Z, so that the result looks pleasing on the two boundaries, without having to greatly increase the step size of vectors or computing needless values outside from , ?

When I say "focus the camera" I mean the following two lines of code:

xlim([-5 5]);
ylim([-5 5]);