function dxdt = siwrODE(t,x,params)

%x(1) = susceptibles s
%x(2) = infecteds i
%x(3) = pathogen in water w
%x(4) = recovereds r

mu = 0;
betaI = params(1);
betaW = params(2);
xi = params(3);
gamma = 0.25;
%k will be params(4);

dxdt = zeros(4,1); %column vector for the state variables
dxdt(1) = mu - betaI*x(1)*x(2) - betaW*x(1)*x(3) - mu*x(1);
dxdt(2) = betaI*x(1)*x(2) + betaW*x(1)*x(3) - gamma*x(2) - mu*x(2);
dxdt(3) = xi*(x(2) - x(3));
dxdt(4) = gamma*x(2) - mu*x(4);
