function SIR()
	% Model-specific parameters
	N = 1000;
	beta = 0.01;
	gamma = 5;
	mu = 0.1;
	
	% Solver-specific parameters
	Tmax = 20;	% How long to solve for
	dt = 1/100;	% Step length
	
	F_S = @(S, I, R) -beta*S.*I + mu*R;
	F_I = @(S, I, R) beta*S.*I - gamma*I;
	F_R = @(S, I, R) gamma*I - mu*R;
	
	[S,I] = meshgrid(linspace(0,N,20));
	U = F_S(S, I, N-S-I);
	V = F_I(S, I, N-S-I);
	U(S+I>N) = 0;
	V(S+I>N) = 0;
% 	quiver(S,I,U,V);
% 	hold on;
% 	plot([0,N], [N,0], 'r');
% 	xlabel('S');	ylabel('I');
% 	return;
	
	
	nt = ceil(Tmax/dt) + 1;
	t = linspace(0, Tmax, nt);
	dt = t(2)-t(1);
	s = zeros(nt,1);
	i = zeros(nt,1);
	r = zeros(nt,1);
	
	% Initial data
	i(1) = 10;
	s(1) = N-i(1);
	r(1) = 0;
	
	% Set up the plotting window
	figure();
	plotInterval = ceil(nt / 50);	% Plot every 'plotInterval' time step
	
	for n = 2:nt
		% Advance the solution one step
		s(n) = s(n-1) + dt*F_S(s(n-1), i(n-1), r(n-1));
		i(n) = i(n-1) + dt*F_I(s(n-1), i(n-1), r(n-1));
		r(n) = r(n-1) + dt*F_R(s(n-1), i(n-1), r(n-1));
		
		% Plot the solution
		if mod(n,plotInterval) == 0 || n==nt
			tiledlayout(2, 2);
			nexttile;
			
			% Plot the phase space
			quiver(S,I,U,V);
			hold on;
			plot(s(1:n), i(1:n), 'r');
			plot(s(n), i(n), 'r*');
			title(sprintf('S = %.1f, I = %.1f', s(n), i(n)));
			xlabel('S');
			ylabel('I');
			hold off;
			
			% Plot each solution up to this time
			plotvars = { s(1:n), i(1:n), r(1:n) };
			varnames = { 'S', 'I', 'R' };
			for k = 1:3
				nexttile;
				plot(t(1:n), plotvars{k});
				xlim([0, Tmax]);
				xlabel('t');
				ylabel(varnames{k});
			end
			drawnow;
		end
	end
end