function rec = GillespieMatlab()
% simple Gillespie simulation of A->B; B->A
% uncommented: exercise is to identify functional aspects of the code

% this vector contains the rates k_1, k_2: modify these later
  c = [2 2];
  
  x = [1 0];
  t = 0; endt = 100; timelast = 0;
  N = 2; M = 2;
  stoichiometry = [-1 1; 1 -1];
  reactants = -stoichiometry.*(stoichiometry < 0);
  a = zeros(M,1);
  rec = zeros(endt,N);
  occ = [0 0];
  
  while(t < endt)
      h = ones(M,1);
      for j=1:M
          for k=1:N
              h(j) = h(j)*CombinatoricConcentration(x(k), reactants(j,k));
          end
      end
      
      az0 = 0;
      for j=1:M
          a(j) = h(j)*c(j);
          az0 = az0 + a(j);
      end
      
      r1 = rand(); r2 = rand();
      tau = (1/az0)*log(1/r1);

      tmp = r2*az0;
      cumsuma = cumsum(a);
      
      reaction = -1;
      for j=1:M
          if cumsuma(j) > tmp && reaction == -1
              reaction = j;
          end
      end
         
      timenow = floor(t+tau);
      if(timenow > timelast)
        for tmpt=timelast:timenow
          rec(tmpt+1,:) = x;
        end
        timelast = timenow;
      end
      
      x = x + stoichiometry(reaction,:);
      t = t + tau;
 
  end
 
      
  end
          
  
