@Scot Gould Although @Rouben Rostamian has clarified that the any Robin boundaries can be solved, it seems that Maple finds sometimes the right answer consistent for the initial condition given and sometimes not. Your first solution was already compatible for the right-hand boundary (at 50 points it is pretty close on almost all the right half of the domain). So I knew to play around with the left-hand one, and used an analogy with the solutions that I was most familar with to figure out that sign change. I usually solve cases that evolve toward steady state, so I added my bias there. But maybe Maple is also consistently solving some classes of boundary conditions; it would help if we know more clearly what assumptions it is making. In this respect, @nm's suggestion to put in the initial solution in afterwards might be counterproductive; maybe Maple is not returning a solution because it knows it can't solve the case compatible with that initial conditions. I suspect very specific assumptions on signs of parameters at the pdsolve stage might be helpful.

I don't think the explanation on the Wikipedia page is that helpful; one needs to think of the relative signs to achieve the right physical outcome. I deliberately kept the fluxes the *same* sign to get to a steady-state, then adjusted the other signs. The general steady state solution is a straight line. If the slope of that is positive, the flux is negative and there will be the same heat flux coming in the right-hand side as going out the left-hand side at steady state. So if the initial fluxes are not equal but of the correct sign, then the other parts of the Robin boundary condition have to allow the slopes to evolve to be the same (assuming you want steady-state of course). In your case the fluxes evolve to be zero.

I didn't really want to use all that opcrap to extract the pieces and did try to do it without a loop first, but I also didn't want to spend to much time optimizing it. Make a procedure that determines whether it is an eigenvalue solution without an analytical solution for the eigenvalues might be worthwhile, but in that case I would probably solve it "by hand". Actually with Maple that is not too difficult. If you call pdsolve with HINT=f(t)*g(x), then extracting the ode for t and x and relating them to the eigenvalue is straightfoward. The second-order ode in x can be solved for {g(x),lambda} and if there is an analytical formula for lambda, Maple will find it; otherwise you have to find the eigenvalues numerically.

I should say, I almost never solve the heat/diffusion equation by this eigenvalue method, because the error functions and x/sqrt(k*t) usually found don't easily come out of it. Instead I find the solution by Laplace transforms. (To get the general solution containing x/sqrt(k*t) you need to use SimilaritySolutions.)