I am interested in generating a 1D non-uniform grid on the interval [0, L] with N points, where a region of width $\sigma$ and centred at $\mu$ is at a higher density and where the transition from low and high densities of grid points occurs over a length $\ell$. This is for a finite-difference code, where particular attention is required to a region with large gradients.In my current implementation I specify low and high grid spacings and interpolate between these using a pair of tanh functions. However in this scheme the total number of points isn't known a priori. I can then iteratively adjust the low and high densities until the total number of points is N but this is quite convoluted in practice. My question then is: can anyone help by describing a scheme for generating such a grid in a simpler way, perhaps something analogous to the way you might generate Gauss-Lobato points x(i) = cos(i PI / N) where i is the grid index. 
Generating a non-uniform grid
finite difference;grid
You achieve this with equidistribution to a mesh density function $\rho(x)$.If you consider $x$ as a continuous map from $\xi \in [0,1]$ into your domain $[0,L]$, then the statement $x$ equidistributes $\rho$ is equivalent to$$\int_0^{x(\xi)} \rho(x') \mathrm{d}x' = \xi \theta\,,$$where $$\theta = \int_0^L \rho(x) \mathrm{d}x\,.$$Not sure if this has a name, but I've been calling it MMPDE0, for moving mesh PDE 0 (rather confusing out of context, but from work in Huang, Ren, Russell (1994))  Another way of stating this, is that you're looking for the $x_i$ for $i=0,\dotsc,N$ that satisfy$$\frac1{\theta} \int_{x_{i-1}}^{x_i}\rho(x')\mathrm{d}x' = \frac1N\,.$$In your case, your mesh density function is piecewise linear, so you can find an explicit expression for $\int \rho(x)\mathrm{d}x$ (which will be piecewise quadratic), which you can solve to find the next $x_i$ from $x_{i-1}$.(Note that this will give $N+1$ points in $[0,L]$ with $x_0 = 0$, $x_N=L$, so adjust for your application as necessary).Here is a possible algorithm you could use to find your $x_i$: Here, $y_j$ are the ordered knots of your piecewise linear mesh density function $x_0 = 0$.For $i$ from 1 to $N-1$Find the highest $j$ such that$$ \int_0^{y_j}\rho(x) \mathrm{d}x < \frac{i \theta}{N},$$(trapezium rule)Find $x_i$ such that$$\begin{aligned} \int_0^{x_i}\rho(x')\mathrm{d}x' &=          \int_0^{y_j}\rho(x')\mathrm{d}x' +          \int_{y_j}^{x_i}\rho(x')\mathrm{d}x'\\        &= \frac{i\theta}{N}      \end{aligned}$$(involves solving a quadratic equation in $x_i$)$x_N = L$It is generally a good idea to smooth your mesh density function before equidistributing to it. It should still be okay to represent it as a piecewise linear mesh density function, just with a lot more knots.