I posted the same question on MATLAB central, and Bill Greene kindly provided an answer.
I briefly report and extend it so that it can be useful to others.
Solution by Bill:
Here is one way to create such a plot. Assume you have the point matrix created by the PDE Toolbox mesher, p, and a solution vector, u. The function below will create a plot of that solution along a line defined by the x and y locations of the two end points. My example is for a solution on a unit square and I want a plot along the line (0,.5) to (1,.5). I want to include 25 points in the plot. As you can see, the real work is being done by the TriScatteredInterp function from core MATLAB.
plotAlongLine(p, u, [0,.5], [1,.5], 25);
function plotAlongLine(p, u, xy1, xy2, numpts)
x = linspace(xy1(1),xy2(1),numpts);
y = linspace(xy1(2),xy2(2),numpts);
F = TriScatteredInterp(p(1,:)', p(2,:)', u);
uxy = F(x,y);
figure; plot(x, uxy); //REM: x is chosen here as a curvilinear coordinate
end
I'd like to further remark that the previous function allows one to plot the solution u or its functions f(u), provided u is defined on mesh nodes (as typically happens for FEM approximated solutions).
If one needs to plot sections of functions defined on mesh centers (e.g. functions of grad u), she may use in advance the function pdeprtni which produces node-valued-functions out of center-valued-functions.