function C=model(xr,xp,p,c,t)
	# calculate concentration matrix as a function of time, initial concentrations and model
	x=xp-xr;
	f= @(c,t) transpose(x)*(p.*prod(repmat(transpose(c),rows(xr),1).^xr,2)); #anonymous function describing the ratelaws
	C=lsode(f,c,t);
endfunction
