function [] = KansaLaplace_2D(N,gridtype)
% KansaLaplace_2D
% Script that performs Kansa collocation for 2D Laplace equation
% N: Number of (interior) collocation points, should be a square (e.g. 25,
% 81, 289, 1089)
% gridtype: type of collocation points
%  'u'=uniform, 'h'=Halton, 'c'=Chebyshev
% Calls on: DistanceMatrix

  % IMQ RBF and its Laplacian
  rbf = @(e,r) 1./sqrt(1+(e*r).^2); ep = 3;
  Lrbf = @(e,r) e^2*((e*r).^2-2)./(1+(e*r).^2).^(5/2);
  % Gaussian RBF and its Laplacian
%  rbf = @(e,r) exp(-(e*r).^2); ep = 3;
%  Lrbf = @(e,r) 4*e^2*exp(-(e*r).^2).*((e*r).^2-1);

  % Exact solution and its Laplacian for test problem
  u = @(x,y) sin(pi*x).*cos(pi*y/2);
  Lu = @(x,y) -1.25*pi^2*sin(pi*x).*cos(pi*y/2);
  neval = 40;
  % Compute (interior) collocation points
  switch gridtype
      case 'u'
          [x1, x2] = ndgrid(linspace(0,1,sqrt(N)));
          dsites = [x1(:) x2(:)];
      case 'h'
          dsites = haltonseq(N,2);
      case 'c'
          [x,y] = meshgrid(0.5+0.5*cos(pi*(0:sqrt(N)-1)/(sqrt(N)-1)));
          dsites = [reshape(x,N,1),reshape(y,N,1)];
      otherwise
          error('Please use u, h, or c data types')
  end          
  intdata = dsites;
  % Additional (equally spaced) boundary collocation points
  sn = sqrt(N); bdylin = linspace(0,1,sn)';
  bdy0 = zeros(sn-1,1); bdy1 = ones(sn-1,1);
  bdydata = [bdylin(1:end-1) bdy0; bdy1 bdylin(1:end-1);...
       flipud(bdylin(2:end)) bdy1; bdy0 flipud(bdylin(2:end))];
  % Create additional boundary centers OUTSIDE the domain
  h = 1/(sn-1); bdylin = (h:h:1-h)';
  bdy0 = -h*ones(sn-2,1); bdy1 = (1+h)*ones(sn-2,1);
  bdyctrs = [-h -h; bdylin bdy0; 1+h -h; bdy1 bdylin;...
      1+h 1+h; flipud(bdylin) bdy1; -h 1+h; bdy0 flipud(bdylin)];
  ctrs = [intdata; bdyctrs];
  % Create neval-by-neval equally spaced evaluation locations
  % in the unit square
  grid = linspace(0,1,neval); [xe,ye] = meshgrid(grid);
  epoints = [xe(:) ye(:)];
  % Compute evaluation matrix
  DM_eval = DistanceMatrix(epoints,ctrs);
  EM = rbf(ep,DM_eval);
  exact = u(epoints(:,1),epoints(:,2));
  % Compute blocks for collocation matrix
  DM_intdata = DistanceMatrix(intdata,ctrs);
  LCM = Lrbf(ep,DM_intdata);
  DM_bdydata = DistanceMatrix(bdydata,ctrs);
  BCM = rbf(ep,DM_bdydata);
  CM = [LCM; BCM];
  % Create right-hand side
  rhs = [Lu(intdata(:,1),intdata(:,2)); ...
         u(bdydata(:,1),bdydata(:,2))];
  % Compute RBF solution
  Pf = EM * (CM\rhs);
  % Compute maximum error on evaluation grid
  maxerr = norm(Pf-exact,inf);
  rms_err = norm(Pf-exact)/neval;
  fprintf('RMS error:     %e\n', rms_err)
  fprintf('Maximum error: %e\n', maxerr)
  % Plot collocation points and centers
  figure
  hold on; plot(intdata(:,1),intdata(:,2),'bo');
  plot(bdydata(:,1),bdydata(:,2),'rx');
  plot(bdyctrs(:,1),bdyctrs(:,2),'gx'); hold off
  % Plot solution surface
  figure
  Pfplot = surf(xe,ye,reshape(Pf,neval,neval),reshape(abs(Pf-exact),neval,neval));
  set(Pfplot,'FaceColor','interp','EdgeColor','none')
  [cmin cmax] = caxis;
  caxis([cmin-.25*maxerr cmax]);
  view([-30,30]); 
  colormap hsv
  vcb = colorbar('vert');
  ylim(vcb,[0 maxerr])
  title('Nonsymmetric RBF solution false colored by maximum error.')
  % Plot error surface
  figure
  errorplot = surf(xe,ye,reshape(abs(Pf-exact),neval,neval));
  set(errorplot,'FaceColor','interp','EdgeColor','none')
  [cmin cmax] = caxis;
  caxis([cmin-.25*maxerr cmax])
  view([-30,30]);
  colormap hsv
  vcb = colorbar('vert');
  ylim(vcb,[0 maxerr])
  title('Maximum error for nonsymmetric RBF solution.')
