ODE45 - want to stop simulation when Steady State occurs

So with inspiration from the post "steady state criterion" I have tried to get my simulation to stop when steady state occurs, using an event with ODE45. However it doesn't stop when steady state occurs, it just moves on to the end time defined, which should change to the time of steady state when the script runs. I have two almost similar parts in the script doing almost the same, the second part works great if I fiddle in the numbers missing from the first part, but the first part keeps on getting a t value of 1000 - I know this value should be 30.8 with the current input. Can't seem to figure out why it doesn't work.

These are my scripts:

Inputvalues:

%%%%%%%%%%%%%%%StartScript.m%%%%%%%%%%%%%%%

global yo j c d g b jm cm gm y y2 tf tf2

j= 0.0031;

c=1.27;

d=0.39;

g=0.8;

b=0.5;

jm=0.0042;

cm=0.9525;

gm=0.7600;

yo =[999 1 0];

%%%%%%%%%%%%%%virttest.m%%%%%%%%%%%%%%%

#First part of the script:

global yo y y2 tf tf2

to = 0;

odefun1 = @(t,y) ligeet(y);

efun1 = @(t,y) odeevent1(y);

options = odeset('Events',efun1);

sol = ode45(odefun1,[to 1000],yo,options);

t = (sol.x)';

Y = (sol.y)';

plot(t,Y(:,1),'b',t,Y(:,2),'m',t,Y(:,3),'g',t,Y(:,1)+Y(:,2)+Y(:,3),'k');

tf = t(end);

#second part:

hold on

yt(1)=Y([end],[1]);

yt(2)=Y([end],[2]);

yt(3)=1;

yt(4)=Y([end],[3]);

odefun2 = @(t2,y2) ligeto(y2);

efun2 = @(t2,y2) odeevent2(y2);

options = odeset('Events',efun2);

sol2 = ode45(odefun2,[t(end) 10000],yt,options);

t2 = sol2.x;

t2 = t2';

Y2 = (sol2.y)';

h =plot(t2,Y2(:,1),'b',t2,Y2(:,2),'m',t2,Y2(:,3),'c',t2,Y2(:,4),'g',t2,Y2(:,1)+Y2(:,2)+Y2(:,3)+Y2(:,4),'k');

xlabel('År')

ylabel('Individer')

title('Endemisk SIR model - introduktion af muteret virus med mindsket helbredsrate')

legend(h([1 2 3 4 5]),'S', 'I', 'I*','R','N');

s=[0:70:max(y)];

plot(t(end),s,'.r');

axis([0 t2(end) 0 max(max(y))+min(max(y))]);

tf2 = t2(end);

%%%%%%%%%%%%Ligeet.m%%%%%%%%%%%%%% function F = ligeet(y)

global j c d g b

F(1)=(b*sum(y))-j*y(1)*y(2)-d*y(1);

F(2) = j*y(1)*y(2)-((d+c+g)*y(2));

F(3)=g*y(2)-d*y(3);

F = F(:);

%%%%%%%%%%%%Ligeto.m%%%%%%%%%%%%%%

function F = ligeto(y)

global j c d g b jm cm gm

F(1)=(b*sum(y))-d*y(1)-y(1)*(j*y(2)+jm*y(3));

F(2) = j*y(1)*y(2)-((d+c+g)*y(2));

F(3)=jm*y(1)*y(3)-(d+cm+gm)*y(3);

F(4)=(g*y(2))+ gm*y(3)-d*y(4);

F = F(:);

%%%%%%%%%%odeevent1.m%%%%%%%%%%

function [x,isterm,dir] = odeevent1(y)

dy = ligeet(y);

x = norm(dy) - 1e-5;

isterm = 1;

dir = 0;

%%%%%%%%%%odeevent2.m%%%%%%%%%%

function [x,isterm,dir] = odeevent2(y)

dy = ligeto(y);

x = norm(dy) - 1e-5;

isterm = 1;

dir = 0;

Answer 1

[Venkatachala Sarma]

Hi,

I initially got an error while running your code. I changed the 'y' to 'sol2.y' in the last but 1 and last but 3 lines. Then it started working. I am able to observe the plot only for 1000 and not for 10000 as should have been the case for the second part of the code.

The first plot gives a response for 4 graphs and the second plot gives a few red dots at the maximum value of the solution 2's output.

Are you getting any error? Because in the above code, y is not defined before using it. MATLAB will treat it as a null ([]). So the line in which s is defined, will give an error.

-Venkatachala Sarma