I want to do stop condition which demand dx=0 in coupled differential equations

I have two equations and I want to stop the integration when those two reach to dx=0.

the equations :

dx/dt=2*x-x^2+0.5*x*y

dy/dt=3*y-y^2+0.5*x*y

my code:

alpha=0.5;

beta=0.5;

r1=2;

r2=3;

s1=1;

s2=1;

t0 = 0;

tfinal = 10;

y0 = [1; 1];

[t,y] = ode23(@Two,[t0 tfinal],y0);

yp(1) = (r1+alpha*y(2))*y(1)-s1*(transpose(y(1))*y(1));

yp(2) = (r2 + beta*y(1))*y(2)-s2*(transpose(y(2))*y(2));

plot(t,y)

grid on

title('Two species')

xlabel('time')

ylabel('Population')

legend('X(t)','Y(t)')

function yp = Two(t,y)

yp = diag([2+0.5*y(2)-1*y(1),3+0.5*y(1)-1*y(2)])*y;

end

Answer 1

[Torsten]

alpha=0.5;
beta=0.5;
r1=2;
r2=3;
s1=1;
s2=1;
t0 = 0;
tfinal = 10;
y0 = [1; 1];   
AnonFun = @(t,y)diag([2+0.5*y(2)-1*y(1),3+0.5*y(1)-1*y(2)])*y;
Opt=odeset('Events',@(t,x)myEvent(t,x,AnonFun));
[t,y,te,ye,ie] = ode23(AnonFun,[t0 tfinal],y0,Opt);
te

te = 9.8627

ye

ye = 1×2

4.6658 5.3343

%yp(1) = (r1+alpha*y(2))*y(1)-s1*(transpose(y(1))*y(1));
%yp(2) = (r2 + beta*y(1))*y(2)-s2*(transpose(y(2))*y(2));
plot(t,y)
grid on
title('Two species')
xlabel('time')
ylabel('Population')
legend('X(t)','Y(t)')

function [value, isterminal, direction] = myEvent(t,x,AnonFun)
value = norm(AnonFun(t,x))-1.0e-2;
isterminal = 1;   % Stop the integration
direction  = -1;
end

[shir hartman]

Thank you!!

Answer 2

[Alan Stevens]

You could do something like the following (doesn't actually stop the integration, but only plots the relevant bit):

alpha=0.5;
beta=0.5;
r1=2;
r2=3;
s1=1;
s2=1;
t0 = 0;
tfinal = 10;
y0 = [1; 1];   
[t,y] = ode23(@Two,[t0 tfinal],y0);
yp(1) = (r1+alpha*y(2))*y(1)-s1*(transpose(y(1))*y(1));
yp(2) = (r2 + beta*y(1))*y(2)-s2*(transpose(y(2))*y(2));
plot(t,y)
grid on
title('Two species')
xlabel('time')
ylabel('Population')
legend('X(t)','Y(t)')

function yp = Two(~,y)
if y(1)>14/3
    y = [NaN; NaN];
end
yp = diag([2+0.5*y(2)-1*y(1),3+0.5*y(1)-1*y(2)])*y;

end

(The value of 14/3 comes from setting your two ode's to zero and solving for the resulting value of x.)

[shir hartman]

Thank you , but if I may - I want it to be more specific.

Is there a way to use odeset ? like that :

Opt=odeset("Events",@(t,x)myEvent(t,x,AnonFun));

[t,x]=ode45(AnonFun,[0,5],1,Opt);

function [value, isterminal, direction] = myEvent(t,x,AnonFun)

value = AnonFun(t,x)-1.0e-4;

isterminal = 1; % Stop the integration

direction = -1;

end

(this is from example of one function ...)