Question: 2 domain PDE

I'm working on solving the basic 1D, transient heat transfer equation in a two layer slab configuration.  One side is insulated, the other has a constant heat flux, and there is not contact resistance between the two layers.
 

My formulation looks like this:
equation1:  diff(T1(z,t),t) = k1/rho1/Cp1*diff(T1(z,t),z$2)
equation2:  diff(T2(z,t),t) = k2/rho2/Cp2*diff(T2(z,t),z$2)
IC1:  T1(z,0) = T0
IC2:  T2(z,0) = T0
BC1:  D[1](T1)(0,t) = 0
BC2:  -k2*D[1](T2)(L2+L1,t) = qs
Match1:  -k1*D[1](T1)(L1,t) = -k1*D[1](T1)(L1,t)
Match2:  T1(L1,t) = T2(L1,t)

T1 is solve from 0 to L1, T2 from L1 to L1+L2.

pdsolve doesn't like it because it doesn't recognize the matching conditions as BC's for each domain.

There's an analytical solution but it's nearly unworkable...

Can this be solved in maple?