How do I sort these 2D-Points to get a proper surf plot of this function

Hello! I am currently trying to plot the function

on the unit disk . I wrote a function to get some evaluationpoints inside and on the boundary of the disk but when I plot the function on it I get very weird results. I already discussed a similar problem in a previous question and a nice person explained to me that the problem was the way my points where sorted in the matrix and therefore the reshaping didn't work as desired. Since I get the same problem in the surf plot again, I am quite sure that this is the same case here.

In the function that I wrote to generate points I first started to lay a grid on a square. After that I sort out every point that is outside of the circle or that is on the boundary of the circle. Lastly I put points on the boundary of the circle. The function takes parameters and N. First two parameters are the bounds for the square . The parameter Nis going to be the squareroot of the number of Points we will end up with. The function gives back Matrices xint, xbdy and x. xint holds the interior points and xbdy the points on the boundary of the circle. x should have both points together but there is still the problem on how I sort these... The rows of these matrices are points in that means e.g. x is gonna be a Matrix.

[~,~,~] = GetCircleGrid(-1,1,30); %Example

function [xint,xbdy,x] = GetCircleGrid(a,b,N)
% Radius of Circle
r = (b-a)/2;
% Centerpoint of Circle
x0 = a+r; 
y0 = a+r;

% First get a Grid on the whole square 
[xsqr1,xsqr2] = meshgrid(linspace(a,b,N),linspace(a,b,N));
xsqr = [xsqr1(:),xsqr2(:)];

% Now get rid of Points outside of Circle
k=1; % counter
for j=1:N^2
    if sqrt((xsqr(j,1)-x0)^2+(xsqr(j,2)-y0)^2) < r % If distance from centerpoint is < r
        xint(k,:) = xsqr(j,:);
        k=k+1;
    end
end

NB = N^2-(k-1); % Number of points for boundary so we end up with N^2 points in total

% Get the Boundarypoints
theta = linspace(0,2*pi,NB);
xbdy = zeros(NB,2);
xbdy(:,1) = x0 + r*cos(theta)';
xbdy(:,2) = y0 + r*sin(theta)';

x = sortrows([xint; xbdy]); % How to sort these...? :(

scatter(xint(:,1),xint(:,2),'ob');
axis equal;
hold on;
scatter(xbdy(:,1),xbdy(:,2),'or');
hold off;
end

The picture show, this is exactly how I want it too look for later usage... Now I wanted to plot the above function (for later usage too...) on the unit disc:

% The function
u = @(x) (5/4)*(1-x(:,1).^2-x(:,2).^2).*sin(pi*x(:,1))+1;

% Get some evaluation points ...
N = 30;
[xint,xbdy,x]=GetCircleGrid(-1,1,N); % ... in the unit disk

% Reshape evaluation points
X = reshape(x(:,1),N,N); % not working well...
Y = reshape(x(:,2),N,N); % not working well...

U = reshape(u(x),N,N); % evaluate the function

% Plot the function as well as the evaluationpoints
figure(1)
surf(X,Y,U,'FaceAlpha',0.5,'EdgeColor','interp');
colormap cool; camlight; lighting gouraud; material dull;
title('f(x,y)')
hold on;
scatter3(X,Y,zeros(N,N),'og');
hold off;

As you can see the plot has some weird lines that are not wanted. How should I sort the points in x to get a proper surf plot? In addition to that I wanted to plot the inner points in a different collor than the boundary points but since the number of points inside and the number of points on the boundary of the disc are not square numbers I don't know how to use reshape for xint and xboundary...

Thank you very much for your help! Kind regard, Max.

Answer 1

[Torsten]

x = linspace(-1,1,1000);
y = linspace(-1,1,1000);
[X,Y] = meshgrid(x,y);
Z = nan(size(X));
idx = X.^2+Y.^2 <= 1;
Z(idx) = 1.25*(1-X(idx).^2-Y(idx).^2).*sin(pi*X(idx))+1;
h = surf(X,Y,Z);
set(h,'edgecolor','none')
colorbar