eqn =
Look carefully at your problem.
n = 2.1;
alpha0 = 0;
syms p alpha
eqn = p^(n+1) - alpha0*p^n + p -(alpha0 + alpha) == 0
So you have sort of, almost kind of, a cubic polynomial, but not really. Raising numbers to non-integrer powers causes problems. There is NO analytical solution. Effectively, we can think of this as a very high order polynomial problem, of degree 31 in fact.) And polynomials higher than order 4 have no solution in algebraic form. This has been known for well over 100 years.
Is there something you could do? If you have the value of alpha, then you could simple use fzero. For example, you might do this:
Psol = @(alpha) fzero(@(p) p^(n+1) - alpha0*p^n + p -(alpha0 + alpha),1)
Psol(3)
So plug in any value of alpha, and fzero will (hopefully) return a solution. You could also use vpasolve, which will survive as long as alpha is given in numerical form. Or fsolve, etc. Another way of looking at this is in the form of a plot. (If we recognize that whenever p is less than zero, all sorts of hell will break lose in terms of complex variables, then we can see that p should be positive. So write the problem by isolating alpha from the rest.
alpha = p + p^3.1
Now we can plot alpha, as a function of p.
fplot(@(p) p + p.^(3.1),[0,3])
xlabel p
ylabel alpha
grid on
Given any value of alpha on the vertical axis, you wish to solve for p. Again, there is NO algebraic solution to this. Yes, if you wanted, you could approximate the inverse relationship. What you want to see is that same relation, but with the axes flipped. And you could easily enough do that. For example...
pint = linspace(0,10);
alphaint = pint + pint.^3.1;
plot(alphaint,pint,'-')
xlabel alpha
ylabel p
grid on
You could now easily enough fit those points with a spline. For example...
pspl = spline(alphaint,pint);
fnval(pspl,5)
So when alpha is 5, p would now directly be predicted as 1.4982. Of course, this is just an interpolation. And again, there is no algebraic function you can write down. Sorry.