i attached the code, i erased many lines so you can focus on the relevant lines.
the problem is in line: dydz(1) = -(gradient(dens_mix,z))*(y(1)/dens_mix);
while the parameter is from: dens_mix= partial_dens(1)+partial_dens(2)+partial_dens(3)+partial_dens(4);
and:
partial_dens = [0 0 0 0 ];
for i = 1:length(partial_dens)
partial_dens(i) = (partial_press(i)*molecular_weights(i))/(Rconstant*y(11));
end
the parameters y(11) and partial_press(i) -values are changing along z axis while the others are constants