As @whuber pointed out, statistical methods do not exactly work here. You need to infer the distribution from other sources. When you know the distribution you have a non-linear equation solving exercise. Denote by $f$ the quantile function of your chosen probability distribution with parameter vector $\theta$. What you have is the following nonlinear system of equations:
\begin{align*}
q_{0.05}&=f(0.05,\theta) \\\\
q_{0.5}&=f(0.5,\theta) \\\\
q_{0.95}&=f(0.95,\theta)\\\\
\end{align*}
where $q$ are your quantiles. You need to solve this system to find $\theta$. Now for practically for any 3-parameter distribution you will find values of parameters satisfying this equation. For 2-parameter and 1-parameter distributions this system is overdetermined, so there are no exact solutions. In this case you can search for a set of parameters which minimizes the discrepancy:
\begin{align*}
(q_{0.05}-f(0.05,\theta))^2+ (q_{0.5}-f(0.5,\theta))^2 + (q_{0.95}-f(0.95,\theta))^2
\end{align*}
Here I chose the quadratic function, but you can chose whatever you want. According to @whuber comments you can assign weights, so that more important quantiles can be fitted more accurately.
For four and more parameters the system is underdetermined, so infinite number of solutions exists.
Here is some sample R code illustrating this approach. For purposes of demonstration I generate the quantiles from Singh-Maddala distribution from VGAM package. This distribution has 3 parameters and is used in income distribution modelling.
q <- qsinmad(c(0.05,0.5,0.95),2,1,4)
plot(x<-seq(0,2,by=0.01), dsinmad(x, 2, 1, 4),type="l")
points(p<-c(0.05, 0.5, 0.95), dsinmad(p, 2, 1, 4))

Now form the function which evaluates the non-linear system of equations:
fn <- function(x,q) q-qsinmad(c(0.05, 0.5, 0.95), x[1], x[2], x[3])
Check whether true values satisfy the equation:
> fn(c(2,1,4),q)
[1] 0 0 0
For solving the non-linear equation system I use function nleqslv
from package nlqeslv.
> sol <- nleqslv(c(2.4,1.5,4.3),fn,q=q)
> sol$x
[1] 2.000000 1.000000 4.000001
As we see we get the exact solution. Now let us try to fit log-normal distribution to these quantiles. For this we will use the optim
function.
> ofn <- function(x,q)sum(abs(q-qlnorm(c(0.05,0.5,0.95),x[1],x[2]))^2)
> osol <- optim(c(1,1),ofn)
> osol$par
[1] -0.905049 0.586334
Now plot the result
plot(x,dlnorm(x,osol$par[1],osol$par[2]),type="l",col=2)
lines(x,dsinmad(x,2,1,4))
points(p,dsinmad(p,2,1,4))

From this we immediately see that the quadratic function is not so good.
Hope this helps.