[W,fs] = solverfun(@odefun,tinterval,f0);

function derivs = odefun(W,fs,param);
    F_H2 = fs(1);
    F_CO = fs(2);
    F_CO2 = fs(3);


  ...


    thiele = D_r * sqrt(k_s * rho_c / De);
    eff = 3 / (thiele ^ 2) * (thiele * coth(thiele) - 1);

  ....


derivs = [dF_H2dW; dF_COdW; dF_CO2dW; dF_H2SdW; dF_SO2dW; dF_N2dW; dF_H2OdW; dF_S6dW; dTdW];