Matlab code

clear all;
clc;

global b c m k1 k2

%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;

k1 = 1.09; 
k2 = 0.74;



%% data 
x0 = 2*10^5;    %%uninfected bacteria   
y0 = 2*10^5;    %%infected bacteria     
v0 = 2*10^5;    %%free phage

t0 = 0.0;
tfi =15;
N = 1000;

%% solving system
ti = t0;
tf = tfi;
[x,y,v] = rk4_phage(x0,y0,v0,ti,tf,k1,k2,b,c,m,N);

%%vector time
time = ti:(tf-ti)/(N-1):tf;

figure
semilogy(time,x,time,y,time,v);
legend('uninfected E.coli','infected E.coli','free phage');
title('phage dynamics');
xlabel('time in hours');
ylabel('concentration');

function replication = a(t)
replication = 0.1121*exp(0.0634*t);

function z = f1(t,x,v,k1,b)
z = a(t)*k1*x-b*v*x;
%dx_dt

function z = f2(t,x,y,v,k2,b)
z = a(t)*k2*y + b*v*x;
%dy_dt

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

function [x,y,v] = rk4_phage(x0,y0,v0,ti,tf,k1,k2,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),v(i),k1,b);
    k12 = dt * f2(t,x(i),y(i),v(i),k2,b);
    k13 = dt * f3(t,x(i),y(i),v(i),b,c,m);
    
    k21 = dt * f1(t+dt/2,x(i)+k11/2,v(i)+k13/2,k1,b);
    k22 = dt * f2(t+dt/2,x(i)+k11/2,y(i)+k12/2,v(i)+k13/2,k2,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,v(i)+k23/2,k1,b);
    k32 = dt * f2(t+dt/2,x(i)+k21/2,y(i)+k22/2,v(i)+k23/2,k2,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,v(i)+k33,k1,b);
    k42 = dt * f2(t+dt,x(i)+k31,y(i)+k32,v(i)+k33,k2,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