yes this allows you to fit it to any order of polynomial hence any order of accuracy, cheers! Ryan

Yes thanks Doug, I did your method and evaluated u_x at the end points using [u(5)-u(5-b)]/b b:=0.0001 which worked well. For some reason when I put in your code Robert it doesn't seem to evalf and it just gives me the integral as an expression, I see you got an answer though! Cheers guys, Ryan

Thats great thanks!

Thank you so much, that works great. p.s. I think you must have a much better PC!

PDE := diff(u(x, t), t) = g*(diff(u(x, t), x))+diff(u(x, t), x, x) where g := piecewise(t <= 1/80, piecewise(x <= 2/3, k[1], 2/3 < x, k[2]), 1/80 < t and t <= 2*(1/80), 0, 2*(1/80) < t and t < 3/80, piecewise(x <= 2/3, k[1], 2/3 < x, k[2]), 3/80 < t and t < 4*(1/80), 0, 4*(1/80) < t and t < 5*(1/80), piecewise(x <= 2/3, k[1], 2/3 < x, k[2]), 5*(1/80) < t and t < 6*(1/80), 0, 6*(1/80) < t and t < 7/80, piecewise(x <= 2/3, k[1], 2/3 < x, k[2]), 7/80 < t and t < 8*(1/80), 0, 8*(1/80) < t and t < 9/80, piecewise(x <= 2/3, k[1], 2/3 < x, k[2]), 9/80 < t and t < 10*(1/80), 0, 10*(1/80) < t and t < 11/80, piecewise(x <= 2/3, k[1], 2/3 < x, k[2]), 11/80 < t and t < 12*(1/80), 0) (6 cycles of a sawtooth potential turning on and off) IBC := {u(0, t) = 0, u(1, t) = 0, u(x, 0) = exp(-(x-2/3)^2/(2*a^2))/((2*Pi)^(1/2)*a)} where a:=1/100 pds := pdsolve(PDE, IBC, numeric, time = t, range = 0 .. 1, timestep = 1/10000, spacestep = 1/10000, optimize = true) R := pds:-value(output = listprocedure)

this looks like a good plan, however when i try the bivariate value(output=listprocedure) its saying "length of output exceeds limit of 1000000". Is this some limitation on my laptop or something thats easily fixed? Ryan