% Usage: [x,c] = makegrid(N, rho)
%
% Given a grid spacing function rho(x) and the number of points N, return
% the N points x(:) and the constant c satisfying
%
%       x(i) - x(i-1) = c * rho(x(i-1))
%       x(1) = -1
%       x(N) + c * rho(x(N)) = 1
%
% For example, x = makegrid(100, @(x) 1) returns a uniform grid.
%
% IMPORTANT: rho(x) need only be defined for x in [-1,1], but you
%            must have rho(x) > 0 everywhere in this interval.

function [x,c] = makegrid(N, rho)
  c0 = 1/rho(-1)/N; % initial guess
  c = fzero(@(c) gridend(N, rho, c) - 1, c0);
  x(1) = -1;
  for i = 2:N
    x(i) = x(i-1) + c * rho(x(i-1));
  end

  function x = gridend(N, rho, c)
    c = max(c, 0); % never step in the negative direction
    x = -1;
    for i = 2:N+1
      % the min(x,1) is to avoid divergences at large x when rho is unbounded
      x = x + c * rho(min(x,1));
    end