clear
%--- The RBF-QR method is used for approximating a function and its 
%--- first and second derivative (in one dimension)
fun1 = @(x) ones(size(x)); 
fun2 = @(x) 25./(25+(x-0.2).^2);
fun3 = @(x) exp(-(x-0.1).^2);
fun4 = @(x) sin(x.^2)-sin(2*x.^2);

funx1 = @(x) zeros(size(x)); 
funx2 = @(x) -50*(x-0.2)./(25+(x-0.2).^2).^2;
funx3 = @(x) -2*(x-0.1).*exp(-(x-0.1).^2);
funx4 = @(x) 2*x.*cos(x.^2)-4*x.*cos(2*x.^2);

funxx1 = @(x) zeros(size(x)); 
funxx2 = @(x) 200*(x-0.2).^2./(25+(x-0.2).^2).^3- 50./(25+(x-0.2).^2).^2;
funxx3 = @(x) -2*exp(-(x-0.1).^2) + 4*(x-0.1).^2.*exp(-(x-0.1).^2);
funxx4 = @(x) 2*cos(x.^2)-4*cos(2*x.^2)-4*x.^2.*sin(x.^2) + 16*x.^2.*sin(2*x.^2);

%--- Select a function
fun = fun3;
funx = funx3;
funxx = funxx3;

%--- Choose problem size and range of epsilon
N =  30; % The number of node points
Ne= 200;  % The number of evaluation points

epvec=logspace(-2,log10(4),50);

%--- Generate node points. Here, we choose a Chebyshev grid, which has
%--- been shown to be beneficial in the small epsilon range.
theta = (pi:-pi/(N-1):0)';
xk = cos(theta);

%--- Uniform evaluation points
xe = linspace(-1,1,Ne)';

%--- Right hand side and exact solution  
f = fun(xk);  
u = fun(xe); 
ux = funx(xe); 
uxx = funxx(xe); 

%--- Distance matrices for the RBF-Direct approach, which will be used
%--- for comparison
[cc,xx]=meshgrid(xk);
rk2 = (xx-cc).^2;
[cc,xx]=meshgrid(xk,xe);
re2 = (xx-cc).^2;

for k=1:length(epvec)
  ep=epvec(k);
  %--- Compute the interpolant using RBF-Direct
  A = exp(-ep^2*rk2);
  B = exp(-ep^2*re2);
  rbfdir(:,k) = B*(A\f);
  %--- Compute the interpolant and its derivatives using RBF-QR
  [A0,Psi] = RBF_QR_diffmat_1D('0',xe,xk,ep);
  [Ax] = RBF_QR_diffmat_1D('x',xe,Psi);
  [Axx] = RBF_QR_diffmat_1D('xx',xe,Psi);
  rbfqr(:,k) = A0*f;
  rbfx(:,k) = Ax*f;
  rbfxx(:,k) = Axx*f; 
  %--- Compute the errors
  errdir(k) = max(abs(rbfdir(:,k)-u)); 
  errqr(k) = max(abs(rbfqr(:,k) -u));
  errx(k) = max(abs(rbfx(:,k)-ux)); 
  errxx(k) = max(abs(rbfxx(:,k)-uxx)); 
end
%--- Display the errors
figure(1),clf
loglog(epvec,errqr,'b',epvec,errdir,'r--',epvec,errx,'k',epvec,errxx,'m')
axis tight
xlabel('\epsilon','FontSize',15)
ylabel('Maximum error','FontSize',15)
legend('RBF-QR','RBF-Direct','1st deriv','2nd deriv')

%--- Find the best solutions in both cases
figure(2),clf
[mi,pos]=min(errqr);

subplot(1,2,1)
H(1)=plot(xe,rbfqr(:,pos)); 
title('RBF-QR: Best solution')
subplot(1,2,2)
H(2)=plot(xe,rbfqr(:,pos)-u); 
title('RBF-QR: Smallest error')

figure(3),clf
[mi,pos]=min(errdir);

subplot(1,2,1)
H(3)=plot(xe,rbfdir(:,pos)); 
title('RBF-Direct: Best solution')
subplot(1,2,2)
H(4)=plot(xe,rbfdir(:,pos)-u); 
title('RBF-Direct: Smallest error')


