Matlab code
global b c m;
%%a replication coefficient of E.coli
%%given by function later
%%b transmission coefficient of phage
b = 1.0000e-006;
%%c replication coefficient of phage
c = 0.706*15;
%%m decay rate of phage
m = 34.8;
 
%% data 1 
x0 = 2*10^5;
y0 = 2*10^5;
v0 = 2*10^5;
t0 = 0.0;
tfi = 15;
N = 500;
 
%% solving system
ti = t0;
tf = tfi;
[x,y,v] = runge_kutta4(x0,y0,v0,ti,tf,b,c,m,N);
%%vector time
time = ti:(tf-ti)/(N-1):tf;
 
semilogy(time,x,time,y,time,v);
legend('uninfected bacteria','infected bacteria','free phage');
title('phage dynamic');
xlabel('time in hours');
ylabel('concentration');
 
figure(2)
plot(time,x,time,y,time,v);
legend('uninfected bacteria','infected bacteria','free phage');
title('phage dynamic');
xlabel('time in hours');
ylabel('concentration');


function z = f1(t,x,y,v,b)
z = a_rep(t)*y-b*v*x;

function z = f2(t,x,y,v,b)
z = a_rep(t)*y + b*v*x;

function z = f3(t,x,y,v,b,c,m)
z = c*y - b*v*x - m*v;

function [x,y,v] = runge_kutta4(x0,y0,v0,ti,tf,b,c,m,N)
 
t = ti;
dt = (tf-ti)/N;
x = zeros(1,N);
y = zeros(1,N);
v = zeros(1,N);
 
x(1) = x0;
y(1) = y0;
v(1) = v0;
for i=1:(N-1)
    k11 = dt * f1(t,x(i),y(i),v(i),b);
    k12 = dt * f2(t,x(i),y(i),v(i),b);
    k13 = dt * f3(t,x(i),y(i),v(i),b,c,m);
    
    k21 = dt * f1(t+dt/2,x(i)+k11/2,y(i)+k12/2,v(i)+k13/2,b);
    k22 = dt * f2(t+dt/2,x(i)+k11/2,y(i)+k12/2,v(i)+k13/2,b);
    k23 = dt * f3(t+dt/2,x(i)+k11/2,y(i)+k12/2,v(i)+k13/2,b,c,m);
    
    k31 = dt * f1(t+dt/2,x(i)+k21/2,y(i)+k22/2,v(i)+k23/2,b);
    k32 = dt * f2(t+dt/2,x(i)+k21/2,y(i)+k22/2,v(i)+k23/2,b);
    k33 = dt * f3(t+dt/2,x(i)+k21/2,y(i)+k22/2,v(i)+k23/2,b,c,m);
    
    k41 = dt * f1(t+dt,x(i)+k31,y(i)+k32,v(i)+k33,b);
    k42 = dt * f2(t+dt,x(i)+k31,y(i)+k32,v(i)+k33,b);
    k43 = dt * f3(t+dt,x(i)+k31,y(i)+k32,v(i)+k33,b,c,m);
    
    t=t+dt;
    x(i+1) = x(i) + 1/6*(k11+2*k21+2*k31+k41);
    y(i+1) = y(i) + 1/6*(k12+2*k22+2*k32+k42);
    v(i+1) = v(i) + 1/6*(k13+2*k23+2*k33+k43);
    
end