K = 1e-5;  % M
kd=10;     % arbitrary, set to alter convergence speed, numerical stability
dt = 1e-7; % arbitrary, alter to change convergence speed, numerical stability
v0 = [0.2; 0.2 ; K]; % M [HA], [A], [X],

%% for dilute conditions
%kd=80;
%dt = 1e-5;
%v0 = [0.0001; 0.0001 ; K]; % [HA], [A], [X]

Nstep = 10000;
t = [0:dt:(Nstep-1)*dt];
tdilute = t(round(Nstep*[0.1:0.1:0.9]));

dv_dt = @(kd,K,v) [ -kd*v(1)+kd/K*v(2)*v(3); kd*v(1)-kd/K*v(2)*v(3)] ;
dv_ = @(v,dt) dv_dt(kd,K,v)*dt;
dv = @(v) dv_(v,dt);
v = v0;
conc = v0;
for ii=2:Nstep
    vstep = dv(v);
    v = v + [vstep(1); vstep(2); vstep(2)];
    conc(:,ii) = v;
    if any(t(ii) == tdilute)
        v = v/1.2;  % <-- change factor to dilute by here
    end
end

figure
subplot(2,1,1)
plot(t,conc(1,:),'r')
hold on
plot(t,conc(2,:),'g--')
plot(t,conc(3,:),'b--')
ylabel('[c]')
legend('[HA]','[A^-]','[H^+]')
subplot(2,1,2)
plot(t,-log10(conc(3,:)),'b--')
ylabel('pH')
xlabel('time (au)')