Thresholds in ODE solvers

Working with models of ecological communities based on diff. equations, I want to define a extinction threshold. If, during the simulation, any species density drops below that it should go directly to zero. I do it with a line in the equations file like:

x(x < threshold)= 0;

But it's not working. As a simple example see this:

function dn= expdecay(t,x)

dn= zeros(1,1);

dn= -0.1*x;

x(x<10)= 0;

When I call: [t,x]= ode45('expdecay',[1 100],[100]); plot(t,x)

I'd expect a straight drop to zero when x= 10, but it doesn't happen. Any idea how to implement this?

Answer 1

[Teja Muppirala]

Make your expdecay look something like this:

dn= zeros(1,1);
dn= -0.1*x;

dead = x < 10;
assignin('caller','dead',dead);
evalin('caller','y(dead) = 0;');

[Teja Muppirala]

This line is not needed: dn= zeros(1,1);

[Francisco de Castro]

Worked beautifully

Answer 2

[Sean de Wolski]

Your function is returning dn. x isn't passed back so the last line doesn't do anything...