title: Using NonNegative setting in Matlab odeset()
body: <p>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?</p>
  
  <p>Here's some code:</p>
  
  <pre><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)];
  </code></pre>
  
  <p>The code for equation(...) is:</p>
  
  <pre><code>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
  </code></pre>