sol = solve(L0 -1 +a*L1^(a-1) -a*A*L0^(a-1)==0, L1 -1 -a*L1^(a-1) +a*A*L0^(a-1)==0,'returnconditions',true);
sol =
struct with fields:
L0: [28×1 sym]
L1: [28×1 sym]
parameters: [1×1 sym]
conditions: [28×1 sym]
>> sol.conditions
ans =
A == 0
A == 0
A == 0
A == 0
A == 0
A == 0
A == 0
A == 0
A == 0
A == 0
A == 0
A == 0
A == 0
A ~= 0
A ~= 0
A ~= 0
A ~= 0
A ~= 0
A ~= 0
A ~= 0
A ~= 0
A ~= 0
A ~= 0
A ~= 0
A ~= 0
A ~= 0
A ~= 0
A ~= 0
So what you have is two sets of solutions. The first 13 solutions are relevant when A == 0. The first is 0, and the next 12 solutions are the cube of the 12 roots of the polynomial z^12 - 5*z^9 + 9*z^6 - 7*z^3 + 46/27 . The 15 solutions after that are the cubes of the 15 roots of the polynomial z^15 - 5*z^12 - 2*A*z^11 + 9*z^9 + 8*A*z^8 + (4*A^2*z^7)/3 - 7*z^6 - 10*A*z^5 - 4*A^2*z^4 - (z^3*(8*A^3 - 46))/27 + 4*A*z^2 + (8*A^2*z)/3 + (16*A^3)/27 .
If your 
is intended to indicate the open set, then A cannot be 0, and you can
is intended to indicate the open set, then A cannot be 0, and you can syms L0 L1 A
assume(A>0 & A<1)
sol = solve(L0 -1 +a*L1^(a-1) -a*A*L0^(a-1)==0, L1 -1 -a*L1^(a-1) +a*A*L0^(a-1)==0)
to get just the 15 solutions
You can convert that polynomial into a roots() call:
roots([1, 0, 0, -5, -2*A, 0, 9, 8*A, (4*A^2)/3, -7, + -10*A, -4*A^2, 46/27 - (8*A^3)/27, 4*A, (8*A^2)/3, (16*A^3)/27])
It appears to me that the 3rd and 6th and 15th root are real-valued. I am not sure at the moment if there are any particular real points for specific A values; I am seeing artifacts due to the way that roots are numbered.
