% plot the solution
% replace the dummy variables of the solution
% in the extension by NaN (then they are not plotted)
INS = find(M==1); OUT = find(M==0);

% convert to double for Matlab 6
INS1 = double(M); 
INS1(OUT) = NaN;  

UU = reshape(U,GRID.nx,GRID.ny).*INS1;

figure('name','solution')
surf(GRID.X,GRID.Y,UU);
xlabel('x')
ylabel('y')
zlabel('u')
title('solution')
hold on

% add all intersection points
% value at the intersection point is -[u],
% it is, J(1:6:end)
P=plot3(IX.coord(:,1),IX.coord(:,2),-J(1:6:end),'r.');
legend('intersection points')
set(P,'markersize',15)

axis equal
box on