Take the periodicity of the diffequation into account. Since t only occurs in Sin or Cos with minimal 10 t, you have peroidicity of 2 Pi/10. NDSolve for high t would be very unreliable; integrate within lower period (But not from the beginning t==0 because of disturbing effects due to initial conditions). Rationalize parameters. sbar = 420; theta = 1; h = 1/2; w2 = 10; nu = 7/10; al[t_] = sbar Sin[2 w2 t]; rho[t_] = 1 + nu Sin[w2 t]; r[x_, t_] = -theta (1 - x) D[rho[t], t] // Simplify; r1[x_, t_] = -theta x (1 - x) al[t] rho[t] (1 - 2 x + h (4 x - 3)) // Simplify; (eqs = Subtract @@@ {D[v[x, t], t] == -al[t] x (1 - x) D[(x + h (1 - 2 x)) v[x, t], x] + (x (1 - x)/(2 rho[t])) D[v[x, t], x, x] + r[x, t] + r1[x, t], v[x, 0] == 0, v[0, t] == 0, v[1, t] == 0} // Together // Numerator // Expand // Simplify) // TableForm vsol = v /. First@NDSolve[Thread[eqs == 0], v, {x, 0, 1}, {t, 0, 3}, MaxSteps -> 10^5, MaxStepSize -> {.0002, .04}, StartingStepSize -> {.0001, .01}] Periodiity in t for all x shown in graphic. Integrating from t=2 Pi/w2 to t=4 Pi/w2 should yield the same as your given range. Manipulate[ Plot[{vsol[x, t], vsol[x, t + 2 Pi/w2], vsol[x, t + 4 Pi/w2]}, {t, 2 Pi/w2, 4 Pi/w2}, PlotPoints -> 100], {{x, 1/2}, 0, 1}] Plot3D[vsol[x, t], {x, 0, 1}, {t, 2 Pi/w2, 4 Pi/w2}, PlotRange -> All, PlotPoints -> 50] nint1 = NIntegrate[2 vsol[x, t], {x, 0, 1}, {t, 2 Pi/w2, 4 Pi/w2}] (* -0.00674974 *) (* Error message: ....NIntegrate obtained -0.00674974 and 6.36785*10^-8 for the integral \ and error estimates *) If you are not satisfied with the estimated error of 6 10^-8, you have to use higher WorkingPrecision. But this takes very long and may crash due to not enough memory. Integrate[ 2 (1 - x) rho[t], {x, 0, 1}, {t, 42 Pi/w2, 44 Pi/w2}]/(2 Pi/w2) + nint1/(2 Pi/w2) (* 0.989257 *)