interp2 meshgrid problem

[Grant Stephens]

Hi

Running version 3.2.4 on Ubuntu 12.04 and have come across something that doesn't make sense to me.

The following program gives the error that the data is not in meshgrid format but all the xx, yy and xxx, yyy data is generated from the meshgrid command. I suspect it is becuase it is a non-linear meshgrid? It there a workaround or it a limitation?

Thank You

Grant

% Set up grids and tensor product Laplacian and solve for u:

  N = 24; [D,x] = cheb(N); y = x;

  [xx,yy] = meshgrid(x(2:N),y(2:N));

  xx = xx(:); yy = yy(:);       % stretch 2D grids to 1D vectors

  f = 10*sin(8*xx.*(yy-1));

  D2 = D^2; D2 = D2(2:N,2:N); I = eye(N-1);

  L = kron(I,D2) + kron(D2,I);                       % Laplacian 

  figure(1), clf, spy(L), drawnow

  tic, u = L\f; toc          % solve problem and watch the clock

% Reshape long 1D results onto 2D grid:

  uu = zeros(N+1,N+1); uu(2:N,2:N) = reshape(u,N-1,N-1);

  [xx,yy] = meshgrid(x,y);

  value = uu(N/4+1,N/4+1);

% Interpolate to finer grid and plot:

  [xxx,yyy] = meshgrid(-1:.04:1,-1:.04:1);

  uuu = interp2(xx,yy,uu,xxx,yyy,'cubic');

  figure(2), mesh(xxx,yyy,uuu), colormap(1e-6*[1 1 1]);

  xlabel x, ylabel y, zlabel u

  text(.4,-.3,-.3,sprintf('u(2^{-1/2},2^{-1/2}) = %14.11f',value))

 function [D,x] = cheb(N)

  if N==0, D=0; x=1; return, end

  x = cos(pi*(0:N)/N)'; 

  c = [2; ones(N-1,1); 2].*(-1).^(0:N)';

  X = repmat(x,1,N+1);

  dX = X-X';                  

  D  = (c*(1./c)')./(dX+(eye(N+1)));      % off-diagonal entries

  D  = D - diag(sum(D'));