function [t, x] = Poisson_timestep(S, M, h, c, t)
% POISSON_TIMESTEP   Simulates the tau-leap (Poisson approximation)
%
% [t, x] = Poisson_timestep(S, M, h, c, tspan)
%
% Inputs:
% S      Stoichiometry matrix (matrix of size n x v, where 
%        n is the number of species and 
%        v is the number of reactions).
% M      Initial state of the system (vector of length n)
% h      Hazard function. Takes two arguments x (vector of length n)
%        and c (vector of length v) as input and
%        outputs a vector of length v.
% c      Stochastic rate constants (vector of length v).
% t      Time span, a vector giving the time discretisation.
%
% Outputs:
% t      Time discretisation (same as t)
% x      Vector of system states
%
% Note that all vectors are column vectors.
% 
% 2019, Henrik Mannerström

x = zeros(numel(M), numel(t));

% Initialize
current_state = M(:);
x(:, 1) = current_state;

for idx = 2:numel(t)
    %% Calculate the reaction rate vector
    % ... #### ....
    
    %% Compute Delta t
    Delta_t = t(idx) - t(idx - 1);
   
    %% Sample Poisson random numbers
    % ... #### ...
    
    %% Update the current state
    % ... #### ...
    
    %% Truncate to zero
    current_state(current_state < 0) = 0;
    
    %% Store computed value
    x(:, idx) = current_state;
end