How can I obtain ODE solution to higher precision?

I am trying to solve a system of differential equations that represent a chemical reaction. 3 of the functions(representing concentrations) need to have values that are very high relative to the changes they undergo. So at given point y = 0.1 and dy = 10^-53, which results in a solution vector that suggests y is constant. I know typically this issue can be worked around by using the digits() and vpa() fucntion, however I am not sure how I can make the ode solver use this level of presicison. Decreasing the Relative Tolerance does not seem to solve the issue. The equations are stiff so I am not sure if solving them symbolically through dsovle is an option.

Here is my code for reference:

n = 10000;
tspan = [0 logspace(-12,-4, n)];

y0 = [0, 0, 0.05, 0.101, 0.1, 0];
options = odeset('RelTol',1e-8,'AbsTol',1e-10);
[T,unper(:,:)] = ode23s(@CombODE, tspan, y0, options);

function dy = CombODE(t, y, flag)
  K = [4.45740519775544e-08;194796303866372;4.71e+15;3.8e+16]
  Km = [8.46547310572481e+47;1.43069785956689e-49;7.41464500062955e-40;1.83447707114187e-47]
    H = y(1);
    O = y(2);
    O2 = y(3);
    H2 = y(4);
    OH = y(5);
    H2O = y(6);
    
    R1 = K(1)*H2;
    R2 = K(2)*O*O;
    R3 = K(3)*O*H;
    R4 = K(4)*H*OH;
    Rm1 = Km(1)*H*H;
    Rm2 = Km(2)*O2;
    Rm3 = Km(3)*OH;
    Rm4 = Km(4)*H2O;
    
    dy = [2*R1-2*Rm1-R3-R4+Rm3+Rm4;
        -2*R2+2*Rm2-R3+Rm3;    
        R2-Rm2;
        -R1+Rm1;
        R3-Rm3-R4+Rm4;
        R4-Rm4;
        ];
end

[Torsten]

K = Kdiss and Km = kmdiss ?

[Vishnu Jayaprakash]

Yes, sorry I shall update the question to be more clear.