Dear Jos,
You are right that the code that I posted uses Mathematica 10 specific commands.
I will have a look at your code. But for the time being, I rewrote the code so that it runs in Mathematica 9 as well. You cannot run the region functions etc. You still need the img2 and img3 images, and then you execute:
sols = Quiet[
NDSolve[{Div[
1./500.*(diffcoeff[798.*x, 654*y])^4*
Grad[u[t, x, y], {x, y}], {x, y}] - D[u[t, x, y], t] +
0.025*u[t, x, y] == 0,
u[0, x, y] == Exp[-1000. ((x - 0.6)^2 + (y - 0.6)^2)],
u[t, 0, y] == 0, u[t, 1, y] == 0, u[t, x, 0] == 0,
u[t, x, 1] == 0}, u[t, x, y], {t, 0, 20}, {x, 0, 1}, {y, 0, 1}]]
Plotting works like this:
ImageCompose[img3, {ContourPlot[
u[t, x, y] /. sols /. t -> 8, {y, 0, 1}, {x, 0, 1},
PlotRange -> {{0, 1}, {0, 1}, {0.01, All}}, PlotPoints -> 100,
Contours -> 200, ContourLines -> False, AspectRatio -> 798./654.,
ColorFunction -> "Temperature"], 0.6}]
Please do refer to the attached notebook as well. Interestingly, I noted a couple of things:
1) The evaluation takes much (!) longer in Mathematica 9. So when you execute the attached notebook you might want to have a cup of tea.
2) I think that I ran the simulation with the same parameters, but the shape of the tumor is slightly different. This might be because of the different integration schemes. So if you want to use this, you need to investigate it. All comments are very welcome.
Here's the frame that corresponds to the one I put up in the original post.

Obviously, the edges are a bit different. The edges correspond to relatively low concentrations however. Also, looking at the frame that describe the progression, everything seems to look a bit different:

Just for the sake of completeness, here is the animation.

I hope that this helps,
Marco
Attachments: