```function [x,param,values] = coarserefine(f,intrv,theta,N,alpha)
% interpolation
%
% f: function defined on interval intrv = [a b].
% intrv: interval [a b] where f is defined.
% theta: threshold interval [thetar thetac] where thetac < thetar.
%        thetar: refinement threshold.
%        thetac: coarsening threshold.
% N: Initial equally-spaced N centers.
% alpha: global multiplier of multiquadric parameters
%
% Example 1 :  f = @(x) abs(x+.04)
%              coarserefine(f,[-1 1],[5e-5 5e-7],11,0.75)
%
% Example 2 :  coarserefine(@(x)myfun(x,5),[-1 1],[5e-5 5e-7],11,0.75)
%              %----------------------%
%              function y = myfun2(x,c)
%              y = 1./(1 + (c*x).^2);
%              %----------------------%
% For reference, see:
% interpolation and collocation problems. submitted to Computers Math. Appl
%
% Tobin A. Driscoll and Alfa R.H. Heryudono    05/14/2006
% MATLAB 7 is recommended.

% Initial points
x = linspace(intrv(1),intrv(2),N)';
N = length(x); dx = diff(x);
epsilon = alpha*min([Inf;1./dx],[1./dx;Inf]);
y = x(1:N-1) + 0.5*dx;

refine = true;
while any(refine)
A = zeros(N); B = zeros(N-1,N);
for j=1:N
A(:,j) = mq(x,x(j),epsilon(j));
B(:,j) = mq(y,x(j),epsilon(j));
end
lambda = A\feval(f,x);
resid = abs(B*lambda - feval(f,y));

refine = resid > theta(1);
vals = sortrows([[x;y(refine)] [feval(f,x);B(refine,:)*lambda]]);
x = vals(:,1); vals(:,1)=[];
dx = diff(x); epsilon = alpha*min([Inf;1./dx],[1./dx;Inf]);

coarsen = resid(1:N-2) < theta(2) & resid(2:N-1) < theta(2);
coarsen = 1+find(coarsen); x(coarsen) = []; vals(coarsen) = []; epsilon(coarsen)=[];
fprintf(' Removing %i centers.\n',length(coarsen))

N = length(x);
y = x(1:N-1) + 0.5*diff(x);
end

if nargout > 1
param = epsilon;
if nargout > 2
values = vals;
end
end```