% HermiteLaplaceMixedBCTref_2D
% Script that performs Hermite collocation for 2D Laplace equation
% Note: Prog 36 in Trefethen (2000), exact solution not provided
% 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);
  L2rbf = @(e,r) 3*e^4*(3*(e*r).^4-24*(e*r).^2+8)./...
                       (1+(e*r).^2).^(9/2);
  %  Laplacian for test problem
  Lu = @(x,y) zeros(size(x));
  % Collocation points and centers
  N = 289;  [collpts, N] = CreatePoints(N, 2, 'u');
  collpts = 2*collpts-1;
  indx = find(collpts(:,1)==-1 | collpts(:,2)==-1 | ...
              collpts(:,1)==1 | collpts(:,2)==1);
  bdypts = collpts(indx,:);  % uniform boundary points
  intpts = collpts;  intctrs = collpts;
  sn = sqrt(N);  h = 1/(sn-1);
  bdyctrs = bdypts;  bdyctrs = (1+2*h)*bdyctrs;
  ctrs = [intctrs; bdyctrs];
  % Create neval-by-neval equally spaced evaluation locations
  % in the unit square
  M = 1681;  epoints = CreatePoints(M,2,'u');  epoints = 2*epoints-1;
  % Compute evaluation matrix
  DM_int = DistanceMatrix(epoints,intctrs);
  LEM = Lrbf(ep,DM_int);
  DM_bdy = DistanceMatrix(epoints,bdyctrs);
  BEM = rbf(ep,DM_bdy);
  EM = [LEM BEM];
  % Compute blocks for collocation matrix
  DM_II = DistanceMatrix(intpts,intctrs);
  LLCM = L2rbf(ep,DM_II);
  DM_IB = DistanceMatrix(intpts,bdyctrs);
  LBCM = Lrbf(ep,DM_IB);
  DM_BI = DistanceMatrix(bdypts,intctrs);
  BLCM = Lrbf(ep,DM_BI);
  DM_BB = DistanceMatrix(bdypts,bdyctrs);
  BBCM = rbf(ep,DM_BB);
  CM = [LLCM LBCM; BLCM BBCM];
  % Create right-hand side
  NI = size(intpts,1);  NB = size(bdypts,1);  rhs = zeros(NI+NB,1);  
  rhs(1:NI) = Lu(intpts(:,1),intpts(:,2));
  indx = find(bdypts(:,1)==1);  % Gamma_2
  rhs(NI+indx) = 0.2*sin(3*pi*bdypts(indx,2));
  indx = find(bdypts(:,1)<=0 & bdypts(:,2)==1);  % Gamma_4
  rhs(NI+indx) = sin(pi*bdypts(indx,1)).^4;
  % Compute RBF solution
  Pf = EM * (CM\rhs);
  % Plot collocation points and centers
  figure
  hold on; plot(intpts(:,1),intpts(:,2),'bo');
  plot(bdypts(:,1),bdypts(:,2),'rx');
  plot(bdyctrs(:,1),bdyctrs(:,2),'gx'); hold off
  figure
  xe = reshape(epoints(:,1),sqrt(M),sqrt(M));
  ye = reshape(epoints(:,2),sqrt(M),sqrt(M));
  surf(xe,ye,reshape(Pf,sqrt(M),sqrt(M)));
  view(-20,45),  axis([-1 1 -1 1 -.2 1]);
  text(0,.8,.5,sprintf('u(0,0) = %12.10f',Pf(841)))