% laying out a 3x3 grid with length 1
d = linspace(0,1,3);  % d= [0    0.5000    1.0000]
[y,x] = meshgrid(d,d);
[xx,yy]=meshgrid(d,d);  % for contour plots
x = [x(:)];
y = [y(:)];
X=[x(:),y(:)];
%--------------------------------------------------------------------------
% Delaunay Triangularization of the domain
Tes = delaunayn(X);
%--------------------------------------------------------------------------
%  Example for Interpolation (linear)
% XI: randomly chosen locations for the interpolation
XI=[.85 .15
    .55 .25
    .05 .8
    .66 .03
     0   0];
[t,P]=tsearchn(X,Tes,XI) % returns the element in which each XI lies and
       % the weight of the location with respect to the 3 vertices of
       % the element
%--------------------------------------------------------------------------
% Random scalar function for the grid
Psi= [1 3.4 5 4 7 8 2.3 5.6 6];
% The interpolation routine for each element of XI
element_cor=Tes(t,:)  % returns the vertices of the element in which XI is 
              % located
psi_element=Psi(element_cor) % the values of Psi at the vertices for the 
              % element containing XI
psi_xi=dot(P,psi_element,2)  % Finally the values of Psi at XI
%--------------------------------------------------------------------------
% Plot for the Grid
%triplot(Tes,x,y)
%--------------------------------------------------------------------------
%  CONTOUR PLOTS --  NOT FINISHED YET
%Combining the interpolated solutions with the solution vector
Psi=[Psi,psi_xi'];
%Psi=reshape(Psi,sqrt(length(Psi)),sqrt(length(Psi)));

% contour(xx,yy,Psi');
%contour(XI,psi_xi);