The "_x

*i*" underscore syntax is only used for equation specification (indicating terms to use partial integration in the FEA discretization) and *can not* be used for postprocessing. Technically it is possible to reconstruct 2nd order derivatives using the

gradient_recovery function as for example here (however accuracy might be an issue, especially at boundaries):

i_dvar = 1; % Dependent variable to differentiate.
j_deriv = 2; % Derivative component 2=x, 3=y, 4=z.
ux = gradient_recovery(fea, i_dvar, j_deriv, fea.sol.u(:,end));
uxx = gradient_recovery(fea, i_dvar, j_deriv, ux);
% Add custom variable for uxx.
fea.vars = struct('name', 'uxx', 'data', uxx);
postplot(fea, 'surfexpr', 'uxx')

Note, that the diffusive flux over the domain does not genrally include second order derivatives as it would be defined as "-d_cd*(cx + cy)" [mol/m^2/s].