I'm interested in your problem.
First of all I guess you know analytical solutions exist for linear heat equations in 1D.
So I suppose your real problem is more complicated than the one you submitted.
For your "toy problem" I began to verify some basic points : for example, "does Maple returns the good solution for still simpler test cases ?".
The simplest one is (1) a single material, with constant thermal conductivity and Dirichlet BCs on the two sides.
Even if the analytical solution is well known, let's just try to verify Maple returns the stationary solution.
You will find in the enclosed file that it is the case : so, without presuming if the transient solution is correctly computed, this is already a good point.
A more complex case is (2) Always a single material, with constant thermal conductivity, but now a Dirichlet BC on the left side and a Robin BC on the right one. I didn't look deeply in the solution Maple gibes, but it seems qualitatively correct.
Now a problem closer to yours (3) Still a single material, with constant thermal conductivity, and Robin Bcs on the two sides (reply to mars sancandi : this problem is well posed ... but the problem with 2 Neuman BCs is not). As you will the solution Maple provides is wrong (look to the argumenty concerning the value of the temperature at the mid point).
This is quite suspect and ominous.
Finally let's try to solve the problem you submitted us.
You write a system with 2 temperatures plus 2 continuity relations at the interface x=L1 (temperature and heat flux continuity).
This is a very common way to proceed when you use Finite Volume of FE methods, because you need these relations to build consistent numerical approximations in the cells (or elements) on the two sides of the alpha discontinuity.
Nevertheless, even with these relations, a good way to proceed is to adapt your mesh accordingly to the thermal diffusion. I do not if you are familiar with this but the idea is to have the same cell Peclet's number in all the cells : simply speaking that means that if alpha = 1 in the left material and 10 in the right one, the length of the cells in the latter has to be 10 times smaller than this same length in the former.
This is somethong pdsolve does not permit.
In the file above you will find the analytical stationary solution for 2 materials with different thermal conduction coefficients and 2 Dirichlet BCs.
This solution is followed by a Maple solution where alpha is defined by a piecewise continuous function. The (quasi) stationary solution is obviously false.
So is the Maple solution obtained with a smooth variant of this previous piecewise continuous alpha, which prooves the discontinuity is not an explanation of the bad result Maple returns.
Finally I end my file by a still simpler problem where alpha is a linear function of x.
Here again the result is obviously false.
Looking to the pdsolve help page one sees that pdsolve can be called in different ways.
It seems to me (but the opinion of someone from the Maplesoft team would be welcome) that it is not possible to use the `numeric`option with a PDE system that contains non partial differential equations (such as your two legitimate continuity constraints).
Unfortunately I have not enough time to spend on you (interesting) problem and I'm not sure Maple can provide the good solution to it. Nevertheless I keep an eye on it and on its developments to come
Hope I helped ?