Plotting derivative of NDSolve solution

I am trying to plot the derivative of an NDSolve solution. My code is :

K = 0.01; u1 = 0; u2 = 1; u3 = 3; xmin = -750; xmax = 750; steep = \
50; tmin = 0; tmax = 500; Diff = 10.;

sol = NDSolve[{D[u[t, x],
t] == -K*(u[t, x] - u1)*(u[t, x] - u2)*(u[t, x] - u3) +
Diff*D[u[t, x], x, x],
u[t, xmin] == u1, u[t, xmax] == u3,
u[0, x] == (u3 + u1*Exp[-x/steep])/(1 + Exp[-x/steep])},
u, {x, xmin, xmax}, {t, tmin, tmax},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MaxPoints" -> 1000, "MinPoints" -> 1000}}]

Animate[Plot[D[u[t, x], x] /. sol, {x, xmin/2, xmax/2},
PlotRange -> {{xmin/2, xmax/2}, {u1 - 0.5, u3 + 0.5}}], {t, tmin,
tmax}]


and I get the errors:

General::ivar: -374.969 is not a valid variable. >>

General::ivar: -374.969 is not a valid variable. >>

General::ivar: -343.719 is not a valid variable. >>

General::stop: Further output of General::ivar will be suppressed during this calculation. >>

General::ivar: -374.969 is not a valid variable. >>

General::ivar: -343.719 is not a valid variable. >>

General::ivar: -312.469 is not a valid variable. >>

General::stop: Further output of General::ivar will be suppressed during this calculation. >>

General::ivar: -374.969 is not a valid variable. >>

General::ivar: -374.969 is not a valid variable. >>

General::ivar: -343.719 is not a valid variable. >>

General::stop: Further output of General::ivar will be suppressed during this calculation. >>

General::ivar: -374.985 is not a valid variable. >>

General::ivar: -374.985 is not a valid variable. >>

General::ivar: -359.679 is not a valid variable. >>

General::stop: Further output of General::ivar will be suppressed during this calculation. >>


I have tried to look at some other answers to this question, all suggesting using the method I have tried.

Need to use Evaluate. Also, not good idea to use UpperCase, K for example is Mathematica symbol. And Manipulate is better to use than Animate. And better to define the pde, ic, bc, separately than put them all inside the call to NDSolve so it is more clear.

k = 0.01;
u1 = 0;
u2 = 1;
u3 = 3;
xmin = -750;
xmax = 750;
steep = 50;
tmin = 0;
tmax = 500;
diff = 10.;

pde = D[u[t, x], t] == -k*(u[t, x] - u1)*(u[t, x] - u2)*(u[t, x] - u3)
+ diff*D[u[t, x], x, x];

bc = {u[t, xmin] == u1, u[t, xmax] == u3};
ic = u[0, x] == (u3 + u1*Exp[-x/steep])/(1 + Exp[-x/steep]);

sol = NDSolve[{pde, ic, bc}, u,{x, xmin, xmax},{t, tmin, tmax},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MaxPoints" -> 1000, "MinPoints" -> 1000}}]

Manipulate[
Plot[Evaluate[D[u[t, x], x] /. sol], {x, xmin/2, xmax/2},
PlotRange -> {{xmin/2, xmax/2}, {0, .05}}],

{{t, tmin, "time"}, tmin, tmax, .1},
TrackedSymbols :> {t}
]