C1=2; C2=1; d1=1; d2=2; b1=0; b2=0;
odefun=@(t,y)[C1*y(1)*(1-b1*y(1)-d2*y(2));
    -C2*y(2)*(1-b2*y(2)-d1*y(1))];
tspan=[0,40];
y0=[2,2];
h=0.01;
Nh=round((tspan(2)-tspan(1))/h);
[tn,un]=rk4(odefun,tspan,y0,Nh);

options=odeset('RelTol',1.e-6);
[tn23,un23]=ode23(odefun,tspan,y0,options);
[tn45,un45]=ode45(odefun,tspan,y0,options);

figure(1)
plot(tn,un(:,1),'b',tn,un(:,2),'b--')
hold on
plot(tn23,un23(:,1),'r',tn23,un23(:,2),'r--')
plot(tn45,un45(:,1),'c',tn45,un45(:,2),'c--')

xlabel('t'); 
legend('RK4 y1','RK4 y','ode23 y1', 'ode23 y2',...
    'ode45 y1', 'ode45 y2')
