You are using a version before R2016b. You will need to code as
DG_dx = bsxfun(@times, (-C-A), (rho-rho_fuel));
The result will be a 100 x 100 matrix, which you will then (:) and return. ode45 will then complain that it is not the same size as the boundary conditions.
Perhaps what you want is instead
DG_dx=(-C-A).'.*(rho-rho_fuel);