A beta distribution seems to suit your needs, but you'll have to perform a transformation in order to change its $(0,1)$ (finite) support to $(-1,1)$ support.
Let $X$ be distributed with a beta distribution, then the random variable $Y$ given by the transformation
$$Y=(b-a)X+a$$
is beta distributed and the PDF has finite support in $(a,b)$. In your case, $a=-1$ and $b=1$. The PDF of this linear transformation is given by:$$p(Y=y|\alpha,\beta,a,b)=f\left(\frac{y-a}{b-a}\right)\frac{1}{b-a},$$
where $f(x)$ is the PDF of the beta distribution given in the wiki page that I cited, and $\alpha$ and $\beta$ are it's parameters. In your case, with $a=-1$ and $b=1$ we get:
$$p(Y=y|\alpha,\beta)=\frac{1}{2}f\left(\frac{y+1}{2}\right).$$
I am guessing that you are looking for a positive, continuous probability distribution with infinite mean and with a maximum density away from zero.
I thought that by analogy with a Gamma distribution ($p(x) \propto x^a \exp(-x) \, dx$), we could try something with a rational (polynomial) rather than an exponential tail. After a little bit of lazy R and Python (sympy) experimentation, I came up with
$$
p(x) = \frac{1}{2\sqrt{3}\cdot \pi/9} \cdot \frac{x}{1+x^3} \, dx
$$
(I initially tried $p(x) \propto x/(1+x^2)$, but its integral diverges.) $\int_0^\infty p(x) \, dx$ is 1, as required, and $\int_0^\infty x p(x) \, dx$ diverges.
I don't know if this distribution has a name/literature associated with it.
The CDF is available in closed form but is pretty horrible ... (see Python code below ...)
$$
3 \sqrt{3} \left(- \frac{\log{\left(X + 1 \right)}}{6 \pi} + \frac{\log{\left(X^{2} - X + 1 \right)}}{12 \pi} + \frac{\sqrt{3} \operatorname{atan}{\left(\frac{2 \sqrt{3} X}{3} - \frac{\sqrt{3}}{3} \right)}}{6 \pi}\right) + \frac{1}{4}
$$
Without actually trying anything, I would guess that distributions of the form $x^a/(1+x^{a+2})$ will generally have these properties (but the computations will get progressively nastier). Someone with more analysis skills could probably prove a bunch of things.
An extremely knowledgeable colleague identified this as almost the same as a "Beta-Type 2 (m=2/3,n=1/3)" distribution (a Beta-Type 2 distribution has a term of the form $(1+x)^n$ in the denominator rather than the $1+x^n$ given above). You might want to use the Beta-Type 2 instead of my version; since you know what it's called you can search for useful code, or literature on its properties (e.g. here or here or McDonald et al 2013), or cite it in a paper: "Beta-Type 2" sounds so much better than "a distribution that some guy on CrossValidated made up".
... the Beta-Type 2 family with its density as
$$
f(x) = \frac{1}{\textrm{Beta}(m,n)} \frac{x^{m-1}}{(1+x)^{m+n}}
$$
over the support $(0,\infty)$
It is evident that if $m$ is chosen to be > 1, the mode will be away from 0. Also, if $n$ is chosen to be $\leq 1$, then the mean will be infinite. This family will produce [an] uncountable number of models with the property you are looking for.
... If you set $Y=X^3$ in your model, then it becomes a Beta-Type 2 $(m=2/3,n=1/3)$ and you can see that this fits the description I have given above.
They also identified the name of @ThomasLumley's contribution:
... it is called power gamma model or exponential-gamma model.
McDonald, James B., Jeff Sorensen, and Patrick A. Turley. “Skewness and Kurtosis Properties of Income Distribution Models.” Review of Income and Wealth 59, no. 2 (2013): 360–74. https://doi.org/10.1111/j.1475-4991.2011.00478.x.
R code:
f <- function(x) 1/(2*sqrt(3)*pi/9)*x/(1+x^3)
integrate(f, 0, Inf) ## 1 with absolute error < 4e-07
curve(f, from=0, to=10)
Python code (because I'm too lazy to integrate):
from sympy import *
x, n, N = symbols('x,n,N')
n=integrate(x/(1+x**3), (x, 0, oo)) ## 2*sqrt(3)*pi/9
integrate(x**2/(1+x**3), (x, 0, oo)) ## infinite mean
cdf = integrate(1/n*x/(1+x**3), (x, 0, X))
print(latex(cdf))
Best Answer
You may be looking for distribution known under the names of generalized normal (version 1), Subbotin distribution, or exponential power distribution. It is parametrized by location $\mu$, scale $\sigma$ and shape $\beta$ with pdf
$$ \frac{\beta}{2\sigma\Gamma(1/\beta)} \exp\left[-\left(\frac{|x-\mu|}{\sigma}\right)^{\beta}\right] $$
as you can notice, for $\beta=1$ it resembles and converges to Laplace distribution, with $\beta=2$ it converges to normal, and when $\beta = \infty$ to uniform distribution.
If you are looking for software that has it implemented, you can check
normalp
library for R (Mineo and Ruggieri, 2005). What is nice about this package is that, among other things, it implements regression with generalized normally distributed errors, i.e. minimizing $L_p$ norm.Mineo, A. M., & Ruggieri, M. (2005). A software tool for the exponential power distribution: The normalp package. Journal of Statistical Software, 12(4), 1-24.