How to use Eventfunction for this ODE set?

The main variables are R and N.

I want to switch between the functions externally or internally.

I know Event function can help but cannot implement it for this case.

EDIT : (working code based on Jan's template)

clc
clear
close all

y0   = [0,0,0,0]; %random initial nos.
t0   = 0;                  %random no.
tEnd = 10;                 %random no.

opt = odeset('RelTol', 1e-3);

y = y0;
t = t0;

State = 1;

while t(end) < tEnd
    yprime = myFunction(t, y, State);
    
    if yprime(2) > yprime(3)
        State = 1;
        disp('state is changed')
    else
        State = 2;
        disp('state is changed')
    end
    
    
    fcn = @(t,y) myFunction(t, y, State);
    
    opt = odeset('Events', @(t, y) myEvents(t,y,State));
    [at, ay] = ode45(fcn, [t(end), tEnd], y(end, :),opt);
   
    t = cat(1, t, at(2:end));
    y = cat(1, y, ay(2:end, :));
end

%solution check
plot(t,y(:,4),'o') %plot results
hold on
z11 =  15*t.^2 - (40/3)*t.^3;
z22 =  50*t.^2 - 1000*t;
plot(t,z11,t, z22) %plot analytical results



figure(2) 
z1 =  30*t - 40*t.^2;
z2 =  100*t - 1000;
plot(t,z1,t, z2) %plot ode set %note the intersection point


function [yprime] = myFunction(t, y, State)
yprime = [5*t;  30*t - 40*t.^2 ; 100*t - 1000 ]; %random funcs.
if State == 1
    yprime(4) = yprime(2);
    disp('state is 1')
else
    yprime(4) = yprime(3);
    disp('state is 2')
end
end


function [val,isterm,dir] = myEvents(t,y,State)

yprime = myFunction(t, y, State);

val      =  yprime(2) - yprime(3);
isterm   =  1;
dir      =  0;

end

Answer 1

[Jan]

y0   = [1,2,3,4];
t0   = 0;
tEnd = 42;

opt = odeset('Tol', 1e-8);

y = y0;
t = t0;
while t(end) < tEnd
    if -y(2) > y(3)
        State = 1;
    else
        State = 2;
    end
    opt = odeset('Events', @(t, y) myEvents(t, y, State));
    fcn = @(t,y) myFunction(t, y, StaTe);
    
    [at, ay] = ode45(fcn, [t(end), tEnd], y(end, :), opt);
   
    t = cat(1, t, at(2:end));
    y = cat(1, y, ay(2:end, :));
end

function myFunction(t, y, State)
if State == 1
    ... implement the first function
else
    ... implement the second function
end
end

function [value,isterminal,direction] = myEvents(t, y, State)
if State == 1
    ... implement the first function
else
    ... implement the second function
end
end

[Angshuman Podder]

@Jan, thanks so much for taking the time to write this up. I understand most of the code except for the Events function. I cannot understand what it is actually doing, so I commented it out. I feel it is telling the solver ode45 to recalculate the step size. I modified the code to test the algorithm is working or not.

Would love your feedback.

P.S. The differential eqns. are dummy eqns.

clc
clear

y0   = [0.1,0.1,0.1,0.21];
t0   = 0;
tEnd = 42;
opt = odeset('RelTol', 1e-3);
y = y0;
t = t0;

State = 1;


while t(end) < tEnd
    yprime = myFunction(t, y, State);
    
    if yprime(2) > yprime(3)
        State = 1;
    else
        State = 2;
    end
    %opt = odeset('Events', @(t, y) myEvents(t, y, State));
    fcn = @(t,y) myFunction(t, y, State);
    
    [at, ay] = ode45(fcn, [t(end), tEnd], y(end, :),opt);
   
    t = cat(1, t, at(2:end));
    y = cat(1, y, ay(2:end, :));
end

plot(t,y)

function [yprime] = myFunction(t, y, State)
yprime = [5*y(1); 4*y(4) + 77*t ; 30*y(4) - 7*t];
if State == 1
    yprime(4) = yprime(2);
    disp('state is 1')
else
    yprime(4) = yprime(3);
    disp('state is 2')
end
end

% 
% function [val,isterm,dir] = myEvents(fun, State)
% if State == 1
%     ... implement the first function
% else
%     ... implement the second function
% val      =  fun(2) - fun(3);
% isterm   =  1;
% dir      =  0;
% end
% end

[Jan]

The event function does not influence the step size, but detects events. This is used to stop the integration and to change the state.

If e.g. the function to be integrated changes, when a component of the trajectory gets 0, this is detected in the event function.