@acer Thank you, this is what I was looking for. This example is more a stand-in for "How do I work with the output from pdsolve in a meaningful way?" which usually involves differentiating, integrating, and etc with solutions. In my case, I just want to verify the conservation law, and numerical differentiation appears to be a necessary evil to work with expressions like int (u_t^2 + c^2u_t^2) dx. The D operator appears to be the way to accomplish this (although Maple probably computes u_t internally in their implementation of the method of lines).

Specifically, I'm plotting something like

plot(evalf(Int(D[2](sol)(x,t), x=a..b, digits=10, epsilon=1e-5, method=_d01ajc)), t=0..t_final)

which is what I mean by the global form of the conservation law - an integral over the entire domain that's supposed to be constant. The local form is a divergence expression. I'm doing this for a 2nd order quasilinear PDE and I want to see the graph of a constant (even to just 2-3 significant digits). If it's not a constant, which I suspect is the case after the wave breaks, then the graph should vary quite a bit. If this is not the case for this conservation law, then perhaps the conservation law int (u_t^2 + c^2u_t^2) dx is violated, by how much, etc.

I also want to plot the global form of the conservation law over time.