x(1)= 600;
y(1)= 60;

gB = 1;
sterberateB = [0.03 0.05 0.075 0.09];
gR = 0.002;
sR = 1;


for sB = sterberateB

f=@(t,x,y) gB*x-sB*x*y;
g=@(t,x,y) gR*x*y-sR*y;

end

schrittweite = [0.2 0.1 0.05 0.001];

for h = schrittweite
t = 0:h:80;

    for i=1:(length(t)-1)
        k1=f(t(i),x(i),y(i));
        l1=g(t(i),x(i),y(i));

        k2=f(t(i)+h/2,(x(i)+0.5*k1*h),(y(i)+(0.5*l1*h)));
        l2=g(t(i)+h/2,(x(i)+0.5*k1*h),(y(i)+(0.5*l1*h)));

        k3=f(t(i)+h/2,(x(i)+0.5*k2*h),(y(i)+(0.5*l2*h)));
        l3=g(t(i)+h/2,(x(i)+0.5*k2*h),(y(i)+(0.5*l2*h)));

        k4=f(t(i)+h, (x(i)+k3*h), (y(i)+l3*h));
        l4=g(t(i)+h, (x(i)+k3*h), (y(i)+l3*h));

        x(i+1)= x(i) + h*(k1 + 2*k2 + 2*k3 +k4)/6;
        y(i+1)= y(i) + h*(l1 + 2*l2 + 2*l3 +l4)/6;


    end
plot(t,x,t,y)
axis([0 100 1 3000])
hold on
end

clear