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);
	tmp_xr=size(xr);
	cfunc= @(t,c) transpose(x)*(p.*prod(repmat(transpose(c),tmp_xr(1),1).^xr,2)); %anonymous function describing the ratelaws
	[junk,C]=ode45(cfunc,t,c);
end

