MATLAB: Switch between 3 different ODE using event functions

eventeventsMATLABodeode45

Hello, I have a system of 3 different equations of motion of second order which I set up as ODEs. These ODEs should toggle on and off depending on the time dependent parameters y(1) & y(3) as seen below, using 1 or more event functions. But this doesn't work with my current code. It only calculates the last of my 3 ODEs. If I delete the last one it calculates only the 2nd one and if I delete the 2nd and 3rd one it doesn't seem to work at all (which is another problem). I need it to automatically switch back and forth between the 3 ODEs. In the end I need a continous plot of my desired variables as also seen below.
Note: I only now that eqmi is the starting equation but I don't know which of the other 2 equation comes next and the time at which the event occurs. It also can switch back and forth between the 3 ODEs multiple times in one timespan.
This was my first question regarding my problem and might help understanding what I'm actually doing, just theoretically, without events. This was my 2nd question, where the answers helped me a lot and shows some of the progress I made but it still doesn't work as needed.
Can anyone help me out with this? Down below is my current Matlab-Code.
function mfc
clc
%%parameters
    m1=10; m2=10; m3=500; k3=300; k2=(k3*(m1+m2))/m3; my=m1/(m1+m2); FW=10; s=0.4; r=0.4; om=5.73; dv0s=2; dv0=dv0s-my*r*om;
%%starting parameters
      tspan = [0 100];
      tstart = tspan(1);
      tend = tspan(end);
      y0 = [0 dv0 0 0];
      tout = tstart;    % copied from ballode example
      yout = y0;
      teout = [];
      yeout = [];
      ieout = [];
   while tout < tend   % main loop
      [t,y,te,ye,ie] = ode45(@(t,y) eqmi(t,y), [tstart tend], y0, odeset('Events', @ZustandI));
      [t,y,te,ye,ie] = ode45(@(t,y) eqmii(t,y), [tstart tend], y0, odeset('Events', @ZustandII));
      [t,y,te,ye,ie] = ode45(@(t,y) eqmiii(t,y), [tstart tend], y0, odeset('Events', @ZustandIII));
      nt = length(t);        % Accumulate output. Copied from ballode example
      tout = [tout; t(2:nt)];
      yout = [yout; y(2:nt,1:4)];
      teout = [teout; te];
      yeout = [yeout; ye];
      ieout = [ieout; ie];
      y0 = [y(nt,1) y(nt,2) y(nt,3) y(nt,4)];  % new IC and tspan
      tstart = t(nt);
  end
%%Plots
      figure(1);                             % displacement
      plot(teout,yeout(:,1),'b--',teout,yeout(:,3),'r:')
      xlabel('time');
      ylabel('displacement');
      title('mfc');
      legend('v2','v3');
      figure(2);                             % velocity
      plot(teout,yeout(:,2),'r:',teout,yeout(:,4),'b--')
      xlabel('time');
      ylabel('velocity');
      title('mfcvelocity');
      legend('v2dot','v3dot');
      figure(3);                             % v2-v3
      plot(teout,yeout(:,1)-yeout(:,3),'r:')
      xlabel('time');
      ylabel('v2-v3');
      title('mfcunterschied');
      legend('v2-v3');
%%functions
   function dy=eqmi(t,y)
      dy=zeros(4,1);
      dy(1)=y(2);                                             % v2dot


      dy(2)=my*r*om^2*sin(om*t);                              % Gl. 2a
      dy(3)=y(4);                                             % v3dot


      dy(4)=(FW*sin(om*t)-k3*y(3))/m3;                        % Gl. 3a
  end
function dy=eqmii(t,y)
    dy=zeros(4,1);
    dy(1)=y(2);                                             % v2dot
    dy(2)=my*r*om^2*sin(om*t)-(k2*(y(1)-y(3)-s))/(m1+m2);   % Gl. 2b II
    dy(3)=y(4);                                             % v3dot
    dy(4)=(FW*sin(om*t)-k3*y(3)+k2*(y(1)-y(3)-s))/m3;       % Gl. 3b II
end
function dy=eqmiii(t,y)
    dy=zeros(4,1);
    dy(1)=y(2);                                             % v2dot
    dy(2)=my*r*om^2*sin(om*t)-(k2*(y(1)-y(3)+s))/(m1+m2);   % Gl. 2b III
    dy(3)=y(4);                                             % v3dot
    dy(4)=(FW*sin(om*t)-k3*y(3)+k2*(y(1)-y(3)+s))/m3;       % Gl. 3b III
end   
%%Events
function [value,isterminal,direction] = ZustandI(t,y)
value = [double(y(1)-y(3)>-s),double(y(1)-y(3)<s)];
isterminal = [1,1];                                     
direction = [-1,1];
end
function [value,isterminal,direction] = ZustandII(t,y)
value = double(y(1)-y(3)>-s);
isterminal = 1;                                     
direction = 0;                                      
end
function [value,isterminal,direction] = ZustandIII(t,y)
value = double(y(1)-y(3)<s);
isterminal = 1;                                     
direction = 0;                                     
end
end

Best Answer

[t,y,te,ye,ie] = ode45(@(t,y) eqmi(t,y), [tstart tend], y0, odeset('Events', @ZustandI));
[t,y,te,ye,ie] = ode45(@(t,y) eqmii(t,y), [tstart tend], y0, odeset('Events', @ZustandII));
[t,y,te,ye,ie] = ode45(@(t,y) eqmiii(t,y), [tstart tend], y0, odeset('Events', @ZustandIII));
Of course this overwrite t,y,te,ye,ie twice. The code computes the complete interval 3 times. But you want to switch between the ODEs.
t   = tstart;
y   = y0;
fcn = @eqmi;
opt = odeset('Events', @ZustandI);
while (t(end) < tend)
  % Run integration until event function stops it:
  % [EDITED, 2018-07-30]: "t" -> "t(end)
  [at, ay, ate, aye, aie] = ode45(fcn, [t(end), tend], y(end, :), opt);
    % Append the new trajectory:
    t = cat(1, t, at(2:end));
    y = cat(1, y, ay(2:end, :));
    % Decide for new equation and event function:
    if y(end, :) > ???? 
      fcn = @eqmi;
      opt = odeset('Events', @ZustandI);
    elseif y(end, :) < ???
      fcn = @eqmii;
      opt = odeset('Events', @ZustandII);
    else
      fcn = @eqmiii;
      opt = odeset('Events', @ZustandIII);
    end
  end
I have no idea how to decide for the specific state, but the idea is to integrate only one each time until the event function stops the integration.
Related Question