Write a function for these coupled ODEs. Inside that function you can use all the normal programmig tools and tricks you can think of. Something like this should work:
function dOISdt = your_chemistry_ode(t,OIS,kP,kD)
O = OIS(1:end-2);
I = OIS(end-1);
S = OIS(end);
dSdt = kD*I;
dIdt = kP*sum(O(1:end-1)) - kD*I;
dOdt = zeros(size(O(:)))
dOdt(1) = -kP*O(1); % Not exactly clear how you want to treat the first spieces in O
for i1 = 2:numel(O)
dOdt(i1) = kP*(O(i1-1)-kP*O(i1));
end
dOISdt = [dOdt;dIdt;dSdt];
end
Then it should be rather straightforward to integrate this with any of the odeNN-functions:
%% Initial conditions
% I just mock something up, completely without ideas about what chain of
% 100 species decaying towards I and S you model
O0 = [100;rand(99,1)];
I0 = pi;
S0 = 0;
OIS0 = [O0;I0;S0];
kP = 1e6; % reaction-rate, no idea about time-constants either
kD = 1e4;
t_span = linspace(0,1,1e6+1); % same again, steps of ms for 1 s...
[tOut,OISout] = ode45(@(t,IOS) your_chemistry_ode(t,OIS,kP,kD),t_span,OIS0);
HTH
