function [ex_sol,ez_sol,A_sol,O_sol,Ab_sol,phi]=solve_biofilm_system(t_max,H,delta_t,delta_z,ex0,A0,O0,Ab0,O_bound,D_o,D_a,eb,delta_ox,delta_xa,delta_xphi,delta_ax,beta_a,t_bar,epsilon,Q,vb,ew,delta_a,delta_ab)

%set up empty matrix's to input the numerical solutions in time and space for each variable
tot_time_steps = t_max/delta_t;
tot_space_steps = (H+vb)/delta_z;
[ex_sol,ez_sol,A_sol,O_sol,Ab_sol] = deal(zeros(tot_space_steps+1, tot_time_steps+1));
phi = zeros(tot_space_steps+1,1);
O0 = zeros(tot_space_steps+1,1);

%Enter values for phi throught biofilm space
for i = 1:(tot_space_steps+1)
    if i <= (H-epsilon)/delta_z
        phi(i) = 0;
    elseif i > H/delta_z + 1
        phi(i) =1;
    else
        phi(i) = 1/epsilon*((i-1)*delta_z-H)+1;
    end
end

%Replace the first column of the matrix for each variable with the initial
%condition
ex_sol(1:H/delta_z+1,1) = ex0;
ez_sol(1:H/delta_z+1,1) = 1-ew-ex0;
A_sol(:,1) = A0;
Ab_sol(:,1) = Ab0;
O = oxygen_solve(H,vb,delta_z,O_bound,D_o,delta_ox,ex_sol(:,1));
O_sol(:,1)=O;
A=A_sol(1:(H/delta_z+1),1);
for j = 1:tot_time_steps
    for i = 1:(tot_space_steps+1)
        ex = ex_sol(i,j);
        ez = ez_sol(i,j);
        O_point = O(i);
        Ab = Ab_sol(i,j);
        %differential equations
        if (i<=H/delta_z+1)
            A_point = A(i);
            ex_sol(i,j+1) = ex + delta_t*(-delta_xa*A_point*(ex) - delta_xphi*phi(i)*(ex));
            ez_sol(i,j+1) = ez + delta_t*(delta_xa*A_point*(ex) + delta_xphi*phi(i)*(ex));
        else
            Ab_sol(i,j+1) = Ab + delta_t*(beta_a*phi(i)*heaviside(t_bar-j*delta_t) - delta_a*Ab + Q/10000*(A_sol(H/delta_z+1,j)-Ab));
        end
        %conditions to ensure that values don't become negative
        if ex_sol(i,j+1)<0
            ex_sol(i,j+1) = 0;
        end
        if A_sol(i,j+1)<0
            A_sol(i,j+1) = 0;
        end
    
    end
    O = oxygen_solve(H,vb,delta_z,O_bound,D_o,delta_ox,ex_sol(:,j+1));
    O_sol(:,j+1)=O;
    A = antimcrobial_solve_release(H,delta_z,D_a,ew,delta_ax,beta_a,Q,ex_sol(:,j+1),Ab_sol(:,j+1),phi,j,delta_t,t_bar,delta_a,delta_ab,eb);
    A_sol(1:(H/delta_z+1),j+1) = A;

end