Rcpp NumericalMatrix data type, global declaration

I'm fairly new with C++ and I'm trying to use it via Rcpp to speed up my R code.

The below code integrates from from t0 to t1- this is done in the "lorenz" function. Test4 integrates using "lorenz" "counts" number of times. However at time "t1" the state of the system is modified in "write_lorenz" before the system is rerun and this is where the problem is. If I run the same program over and over again by calling test4 from R, printing to the screen always produces the same result, however, my returned matrix "u" does not, and seems to eventually converge to whatever "t1" is which is the problem.

My input values don't change so I'm wondering if there a memory leak, or if something else is going on, how to fix it. Also I'm wondering if my initialization of "u" is incorrect and I should be using the "new" command.

What I tried was NumericMatrix* u = NULL; *u = new NumericMatrix; and then I tried accessing the elements of the matrix as *u(1,2) for instance, but accessing the elements this way caused an error saying u was not a function.

Any help would be greatly appreciated

I modified this code from the following site

http://ift.tt/1ChcgjU

so I could use it with Rcpp

//###################################R Code ###############################

library(Rcpp)
sourceCpp("test4.cpp")

sigma  <- 10.0
R <-28.0
b <- 8.0 / 3.0


a2 <-  c(10.0 , 1.0 , 1.0) #initial conditions  X0,I0,Y0
L2 <- c(0.0 , 2.0 , 0.1)  #initial time, kick time, error
counts <- 2 
kick <-1.0; # kick size

pars <-c(sigma,R,b,kick)


test4(a=a,L2=L2,counts=counts,pars= pars) 




// C ++ code 

//[[Rcpp::depends(BH)]]
//[[Rcpp::depends(RcppEigen)]]
//[[Rcpp::plugins("cpp11")]]

#include <Rcpp.h>
#include <RcppEigen.h>
#include <math.h>
#include <boost/array.hpp>
#include <boost/numeric/odeint.hpp>
#include <boost/algorithm/string.hpp>

using namespace boost::numeric::odeint;
using namespace std;
using namespace Rcpp;
using namespace Eigen;


double sigma =0;
double e =0; 
double b =0 ; 
double t0 =0;
double t1 = 0;
double kick =0;

double X0 = 0;
double I0 =0;
double Y0 =0;
double N = 0;

int counter =0;
typedef boost::array< double , 3 > state_type;

NumericMatrix u(4,5);


void lorenz( const state_type &x , state_type &dxdt , double t )
{
    dxdt[0] = sigma * ( x[1] - x[0] );
    dxdt[1] = e * x[0] - x[1] - x[0] * x[2];
    dxdt[2] = -b * x[2] + x[0] * x[1];
}


void write_lorenz( const state_type &x , const double t )
{

  if(t==t1){
  X0 = x[0];
  I0 = x[1];
  Y0 = x[2];
  N = X0+I0+Y0;
  X0 = X0 + exp(kick*N);

  counter++;

//for some reason cout and u don't match for multiple runs of the 
//program
 cout << t << '\t' << X0 << '\t' << I0 << '\t' << Y0 << endl;

  u(counter,0) = t;
  u(counter,1) = X0;
  u(counter,2) = I0;
  u(counter,3) = Y0;




  }

}


// [[Rcpp::export]]
NumericMatrix test4(NumericVector a, NumericVector L2,int counts,NumericVector pars) 
{


double error; // control integration error

// initialize model parameters     
//maybe these should be parenthesis?
sigma = pars[0]; //# average per capita birh rate 1-10(the mean is 4.7)
e = pars[1]; // # per capita average growth rate
b= pars[2];
kick = pars[3]; // kick size

      //cout << sigma << R << b << kick << endl;

     //myfile.open (ST);

  X0 = a[0]; I0 =a[1];  Y0 = a[2];   
  int i = 0;
  t0 = L2[0];
  t1 = L2[1];
  error = L2[2];

u(0,0) =  t0;
u(0,1) = X0;
u(0,2) = I0;
u(0,3) = Y0;


      // initial conditions

  for(i=0;i<counts;i++)
    {

   state_type x = { X0 , I0 , Y0 };
    integrate( lorenz , x , t0 , t1 , error , write_lorenz );
    }


  return u;   // the variable I hope will be global 



}