If you're having trouble getting results closer to the singularity, that's because the value of z at the corner is NaN. Z is generally very large in that vicinity, but it's undefined (and not drawn) at exactly 0,0. The closer you can get the neighboring points, the more of the peak gets represented.
z = @(x,y)(x./(x.^2+y.^2));
x = linspace(0,0.1,100); % just use a finer mesh
y = linspace(0,0.1,100);
[X, Y] = meshgrid(x, y);
surf(X, Y, z(X,Y));
view(-16,18)
colormap(parula)
shading flat
xlabel('x');
ylabel('y');
zlabel('z');
Same would work for a contour map
clf;
contourf(X, Y, z(X,Y),[0 10 20 30 40]); % specific levels
colormap(parula)
shading flat
xlabel('x');
ylabel('y');
One way to avoid needing to use such a large number of points to represent features near the edge or corner of a selected domain is to use logspace instead of linspace. Both these examples use the same number of points.
clf;
% linear mesh spacing
subplot(1,2,1)
x = linspace(0,0.5,20);
y = linspace(0,0.5,20);
[X, Y] = meshgrid(x, y);
h = pcolor(X, Y, z(X,Y));
h.EdgeAlpha = 0.5;
caxis([0 81])
colormap(parula)
xlabel('x');
ylabel('y');
% log mesh spacing
subplot(1,2,2)
n = 20;
x = logspace(0,1.0414,n)/n-1/n;
y = logspace(0,1.0414,n)/n-1/n;
[X, Y] = meshgrid(x, y);
h = pcolor(X, Y, z(X,Y));
h.EdgeAlpha = 0.5;
caxis([0 81])
colormap(parula)
xlabel('x');
ylabel('y');
