How to fix: Index exceeds array bounds?

Hey guys, I have a problem with my code. I programmed an Event-ODE-Function which works fine. Check my previous question fore more background about this here. I post my full code after the jump. Now I want to add another ODE into my existing code which should integrate only once but should change based on the old events and their triggers.

My previous code works alright and I'll mark the new parts with this: %new%. Now I get the error message Index exceeds array bounds. Can you help me fixing it?

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];

    % initiate starting function
    t = tstart;
    y = y0;
    fcn = @eqmi;
    Ifcn = @impulsi;                                      %new%
    opt = odeset('Events', @ZustandI);

    ateout = [];
    ayeout = [];
    aieout = [];
%%main loop
while t(end) < tend
    % Integriere bis zum Stopp durch Event-Funktion
    [at,ay,ate,aye,aie] = ode45(fcn, [t(end) tend], y(end,:), opt);
    [at,I,ate,aIe,aie] = ode45(Ifcn, [t(end) tend], 0, opt);                  %new%

    % new starting parameters
    t = cat(1, t, at(2:end));
    y = cat(1, y, ay(2:end,:));

    ateout = [ateout; ate];
    ayeout = [ayeout; aye];
    aieout = [aieout; aie];

    % new functions and events
    if (-s<y(end,1)-y(end,3) && y(end,1)-y(end,3)<s)
      fcn = @eqmi;
      Ifcn = @impulsi;                                             %new%
      opt = odeset('Events', @ZustandI);
    elseif y(end,1)-y(end,3)<=-s
      fcn = @eqmii;
      Ifcn = @impulsii;                                             %new%
      opt = odeset('Events', @ZustandII);
    elseif y(end,1)-y(end,3)>=s
      fcn = @eqmiii;
      Ifcn = @impulsiii;                                              %new%
      opt = odeset('Events', @ZustandIII);
    end  
end
%%functions
function dy=eqmi(t,y)
    dy=zeros(4,1);
    dy(1)=y(2);                                         
    dy(2)=my*r*om^2*sin(om*t);                             
    dy(3)=y(4);                                          
    dy(4)=(FW*sin(w3*t)-k3*y(3))/m3;                     
end
function dy=eqmii(t,y)
    dy=zeros(4,1);
    dy(1)=y(2);                                           
    dy(2)=my*r*om^2*sin(om*t)-(k2*(y(1)-y(3)-s))/(m1+m2);  
    dy(3)=y(4);                                            
    dy(4)=(FW*sin(w3*t)-k3*y(3)+k2*(y(1)-y(3)-s))/m3;      
end
function dy=eqmiii(t,y)
    dy=zeros(4,1);
    dy(1)=y(2);                                           
    dy(2)=my*r*om^2*sin(om*t)-(k2*(y(1)-y(3)+s))/(m1+m2);   
    dy(3)=y(4);                                             
    dy(4)=(FW*sin(w3*t)-k3*y(3)+k2*(y(1)-y(3)+s))/m3;      
end
function dI=impulsii (t,I)                                  %new%
    dI=k2*(y(1)-y(3)-s);
end
function dI=impulsiii (t,I)                                 %new%
    dI=k2*(y(1)-y(3)+s);
end
function dI=impulsi (t,I)                                   %new%
    dI=0;                                                   % I and dI should be zero at event ZustandI
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 = [0,0];
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

I'll leave the plots out for now. Thanks in advance guys.

[Stephen23]

"Now I get the error message Index exceeds array bounds."

Sure. But you post one hundred lines of code and nothing telling us where the error actually occurs. If you really want to get help, then post the complete error message. This means all of the red text.

[madhan ravi]

W3 undefined

Answer 1

[Torsten]

[at,I,ate,aIe,aie] = ode45(Ifcn, [t(end) tend], 0, opt)

implies that you have only one ODE to solve, but in ZustandI, you address y(3) (which does not exist). That's why you get the error message.

[madhan ravi]

+1 @Torsten scrutinised precisely .

[Lennart]

Thanks. Your help is appreciated as always. I think I get the problem, so I deleted everything in my code above marked : %new% and tried to include all my calculations in the same function, like this:

function dy=eqmi(t,y)
    dy=zeros(6,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(w3*t)-k3*y(3))/m3;                        % Gl. 3a
    dy(5)=y(6);                                             % don't know if I need this at all
    dy(6)=0;                                                % dI
end

and this (+s for eqmiii):

function dy=eqmii(t,y)
    dy=zeros(6,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(w3*t)-k3*y(3)+k2*(y(1)-y(3)-s))/m3;       % Gl. 3b II
    dy(5)=y(6);                                             % don't know if I need this
    dy(6)=k2*(y(1)-y(3)-s);                                 % dI
end

and the starting parameters:

y0 = [0 dv0 0 0 0 0];

but this changes my plots for y(1) to y(4). This cannot be right because dy(5) & dy(6) should be closed calculations in and of it all which cannot affect the other ones. Also it starts right but after plotting like 20 s it changes the plot:

figure(1);                             % displacement
    plot(t,y(:,1),t,y(:,3))
    hold on
    plot(ateout,ayeout(:,1),'rx',ateout,ayeout(:,3),'b.')
    xlabel('time');
    ylabel('displacement');
    title('mfc');
    legend('v2');