I think you are making a silly mistake by considering solution as a List when you have defined it as a InterpolatingFunction object. Here is little modification which may help
Clear[Eq0, FilmEqn, h, Bo, ϵ, K1, δ, Bi, m, r]
Eq0[h_, {Bo_, ϵ_, K1_, δ_, Bi_, m_, r_}] := \!\(
\*SubscriptBox[\(∂\), \(t\)]h\) +
Div[-h^3 Bo Grad[h] +
h^3 Grad[Laplacian[h]] + (δ h^3)/(Bi h + K1)^3 Grad[h] +
m (h/(K1 + Bi h))^2 Grad[h]] + ϵ/(
Bi h + K1) + (r) D[D[(h^2/(K1 + Bi h)), x] h^3, x];
SetCoordinates[Cartesian[x, y, z]];
FilmEqn[Bo_, ϵ_, K1_, δ_, Bi_, m_, r_, time_] :=
Eq0[solution[x, y, time], {Bo, ϵ, K1, δ, Bi, m, r}]
expr = FilmEqn[0, 10^-6, 1, 10^-3, 1, 0.05, 0, time];
fun[a_, b_, t_] := Evaluate[expr /. x -> a /. y -> b /. time -> t];
Plot3D[fun[x, y, 100000], {x, 0, L}, {y, 0, L}, PlotPoints -> 40,Mesh ->None,
ColorFunction -> Hue, PlotRange -> All]
Plot3D[(fun[x, y, #] - fun[x, y, # - 1]), {x, 0, L}, {y, 0, L},
PerformanceGoal -> "Quality",
Mesh -> None, ColorFunction -> Hue,
PlotLabel -> "eq(t)-eq(t-1) at t:= " <> ToString[#]] &
