function R = myfun2 (z, y)
emit = 6*10^-5;
kb = 0.0015;
kh_z = kh_z1 (z);
R = [y(2) ; -kh_z*y(1)+ kb/y(1) + emit^2/y(1).^3];
end




