I have found that MATLAB solves 1D parabolic-elliptic system incorrectly by using pdepe function when Robin boundary conditions are specified.
See, please: https://math.stackexchange.com/questions/829402/matlab-incorrect-computation-with-pdepe-for-parabolic-elliptic-system
If we use Dirichlet boundary conditions, the solution is correct.
Is my reasoning correct?
Some information about this in Russian: http://glebgrenkin.blogspot.ru/2014/05/matlab.html