Saving the output of a function within the pdepe function into an array

I'm trying to save output 's' from a function that's called on by the pdepe partial differential equation solver.

sol = pdepe(m,@pdex14pde,@pdex1ic,@pdex1bc,x,t,options);

Where pdex1ic is an intial conditions vector, pdex1bc is a boundary conditions vector, x is a spatial mesh vector, and t is a time mesh vector.

function [c,f,s] = pdex14pde(x,t,u,DuDx)

Where 's' is defined as

%Reaction:----------------------------------------------- 
s = [(BC0_O2-u(1))*alfax + (-2*R6-0.25*R8) + (-R1)*F;...                    
    R8/F + (-4*R3-2*R7);...                                                 
    -4*R4;...                                                               
    (BC0_SO4-u(4))*alfax + R6 + (-0.5*R5)*F;...                    
    (BC0_Fe-u(5))*alfax - R8 + (4*R3+4*R4+2*R7+R9-R_9)*F;...  
    (BC0_H2S-u(6))*alfax + (-R6-R12) + (0.5*R5-R7+R9-R_9)*F;... 
    -R12/F + (-R9+R_9);... 
    -Ra;...
    -Rb];

Where u(i) are the various chemical species involved in the equations.

'sol' is a three dimensional array of [time x space x chemical species] made up of chemical concentrations at each of those points. These chemical concentrations are a product of the equation c(x,t,u,Du/Dx) * Du/Dt = x^(-m) * D(x^m * f(x,t,u,Du/Dx))/Dx + s(x,t,u,Du/Dx)

I want to generate an array of the same size as 'sol' but only populated with the 's' term values. The time steps taken by the ode15s solver are quite small at parts and don't align with the mesh points for x or t.