You can use the MatLab ODE event location.
I will give an example below, but I will have to make up names since you gave too few info.
Let's assume your ODE function is called 'myfunc' and that E is the third variable
function out = myfunc(t,y)
R = y(1);
Rp = y(2);
E = y(3);
etc.
dR/dt = -k1*R*E+k_1*RE-p1u*CD*R+pd1u*Rp-k3*E7*R+k_3*RV;
dRp/dt = p1u*CD*R-p2u*CE*Rp-k2*Rp*E+k_2*RpE-pd1u*Rp+pd2u*Rpp-k4*E7*Rp+k_4*RpV;
dE/dt = -k1*R*E+k_1*RE-k2*Rp*E+k_2*RpE+p2c*CE*RpE+pi2c*CE*RpEV+k_8*RpEV+k_7*REV;
etc.
out = [ dR/dt;...
dRp/dt;...
dE/dt;...
etc.
]
end
if Ef is the final value of E, your ODE event locator function should look something like this.
function [position,isterminal,direction] = eventFcn(t,y,Ef)
position = y(3)-0.9*Ef; % The value that we want to be zero
isterminal = 1; % Halt integration
direction = 0; % The zero can be approached from either direction
end
You can pass this function to the ODE options structure, which is then pass to the ODE solver
opts = odeset('Events',@(t,y)eventFcn(t,y,Ef));
[t,y] = ode45(@myfunc,tspan,y0,opts);
This should do it for you.
