function SIREquilibrium()
S0=0.8;
I0=0.2;
R=1-S0-I0;
beta=1;
gamma=1/10;
mu=5e-4;

f=@(t,y)[mu-beta*y(1)*y(2)-mu*y(1);beta*y(1)*y(2)-gamma*y(2)-mu*y(2);gamma*y(2)-mu*y(3)];

y1=linspace(0,1,20);
y2=linspace(0,1,20);
[x,y]=meshgrid(y1,y2);
u=zeros(size(x));
v=zeros(size(y));
t=0;
for i=1:numel(x)
    Yprime=f(t,[x(i);y(i)]);
    u(i)=Yprime(1);
    v(i)=Yprime(2);
end
quiver(x,y,u,v,'r')

index exceeds matrix dimensions.

Error in SIREquilibrium>@(t,y)[mu-beta*y(1)*y(2)-mu*y(1);beta*y(1)*y(2)-gamma*y(2)-mu*y(2);gamma*y(2)-mu*y(3)] (line 16)
f=@(t,y)[mu-beta*y(1)*y(2)-mu*y(1);beta*y(1)*y(2)-gamma*y(2)-mu*y(2);gamma*y(2)-mu*y(3)];

Error in SIREquilibrium (line 25)
    Yprime=f(t,[x(i);y(i)]);