Using NonNegative setting in Matlab odeset()

Question

I'm trying to solve some ODEs in MatLab and seeing as the variables in the equations are populations they need to be constrained to being positive. So I tried using odeset() before calling the equation solver to make them non-negative but on plotting the values afterwards they are actually negative at times (in the code below it is the magenta line). What am I doing wrong?

Here's some code:

%Lots of variables
includeJ=1;
cullLIRate=1/2000;
cullDIRate=1/2000;
N = 16800;
beta = 2e-7;
delta = 0.5;
gamma = 1/50;
sigma = 1/400;
mu = 1/365;
maxTime = 30*365;
kappa = N;
gR = 0.05;
mJ = 1/3650;
initJPerAdult = 10;
numInitE = 1000;
TSpan = [0,maxTime];

initState = [N-numInitE,numInitE,0,0,0,initJPerAdult*N];

%IMPORTANT BIT HERE
options = odeset('NonNegative', 1:6)
scirSoln = ode45(@equation,TSpan,initState,[],beta,delta,gamma,sigma,mu,kappa,gR,mJ,cullLIRate,cullDIRate,includeJ);

scirVals = deval(scirSoln,timeToPlot);
plot(timeToPlot,scirVals(1,:));
hold on;
plot(timeToPlot,scirVals(3,:),'k');
plot(timeToPlot,scirVals(4,:),'g');
plot(timeToPlot,scirVals(6,:),'m');

timeToPlot = [0:max(TSpan)/1000:max(TSpan)];

The code for equation(...) is:

function retVal = equation(t,y,beta,delta,gamma,sigma,mu,kappa,gR,mJ,cullLIRate,cullDIRate,includeJ)
retVal = zeros(6,1);

S = y(1);
E = y(2);
LI = y(3);
DI = y(4);
R = y(5);
J = y(6);

retVal(1)= mJ * J - beta * S * (delta * LI + DI);
retVal(2) = beta * S * (delta * LI + DI) - gamma * E;
retVal(3) = gamma * E - (cullLIRate + sigma) * LI;
retVal(4) = sigma * LI - (mu + cullDIRate) * DI;
retVal(5) = mu * DI + cullLIRate* LI + cullDIRate * DI;
retVal(6) = gR * S * (1 - S / kappa) - mJ * J;

end

Answer 1

You are not passing your defined odeset (options variable) to the ODE45 - solver.

The syntax for the ODE45 is: [T,Y] = ODE45(ODEFUN,TSPAN,Y0,OPTIONS,P1,P2...)

Glad it worked! :)