#include <stdlib.h>
#include <stdio.h>
#include "math.h"
#include <mex.h>

#define n      20.0f
#define n_     3.0f
#define mu     0.10f
#define mu_    14.0f
#define gamma  0.10f
#define sigma  0.01f
#define beta   0.50f
#define beta_  0.50f
#define STEPS  0.01f
#define TEND   365.0f
#define COUNT  (int)(TEND / STEPS)
int i;
double yvec[COUNT][2];
void getResults()
{
    yvec[0][0] = 0.0f; yvec[0][1] = 0.0f;
    #define deltay  (double)(sigma*(mu*n + gamma*yvec[i][0]) + beta*(n-yvec[i][0])*(yvec[i][1]/n_) - (mu + gamma)*yvec[i][0])
    #define deltay_ (double)(beta_*yvec[i][0]*(1.0f - (yvec[i][1]/n_)) - mu_*yvec[i][1])
    for (i = 1;i < COUNT; i++) // Euler's Method
    {
        yvec[i+1][0] = yvec[i][0] + STEPS * deltay;
        yvec[i+1][1] = yvec[i][1] + STEPS * deltay_;
    }
    return;
}
double * myvec;
void mexFunction(int nlhs, mxArray *plhs[], 
    int nrhs, const mxArray *prhs[])
{
    getResults();
    plhs[0] = mxCreateDoubleMatrix(COUNT, 2, mxREAL);
    myvec = mxGetPr(plhs[0]);
    myvec = yvec; //this line is wrong, what should it really be?
    return;
}