% model is similar to model 12
% H(phi) in chemotaxis is changed to 0.5*(tanh(phi)+1)
function [N_sol,M_sol,D_sol,G_sol,F_sol,C_sol,S_sol,W_sol,phi_sol,H_sol,A_sol]=solve_metal_ion_wound_healing_bacteria_system(t_max,x_1,x_2,delta_t,delta_x,N0,M0,D0,G0,F0,C0,S0,W0,phi0,H0,A0,N_bound,M_bound,D_bound,G_bound,F_bound,C_bound,S_bound,W_bound,phi_bound,H_bound,A_bound,chi_n,sigma_n,gamma_ncs,kappa_n,gamma_nc,delta_nc,delta_mf,delta_ms,delta_m,beta_gs,delta_gc,beta_gf,S_bar,beta_sc,D_n,D_f,chi_f,sigma_f,rho_b,beta_fb,delta_sn,delta_sf,beta_sm,delta_s,eta,beta_fc,beta_ssf,kappa_f,beta_fcs,gamma_fcs,delta_fc,rho_f,D_c,beta_cb,delta_cn,delta_cf,D_s,beta_sn,beta_ssn,beta_sf,C_bar,gamma_gf,alpha_gmw,alpha_dmw,alpha_fmw,alpha_cmw,beta_wc,beta_wm,beta_wn,beta_sw,delta_w,delta_wf,gamma_gw,N_bar, beta_mn,beta_mg,beta_md,delta_nw,beta_nh,beta_fh,beta_sfh,D_h,delta_hf,delta_hn,delta_h,D_a,delta_aw,delta_a,delta_wa,gamma_b,beta_gfh,sigma_nphi,sigma_fphi,delta_phi,alpha,H_d,H_w,vol_frac,n)

%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 = (x_2-x_1)/delta_x;
[N_sol,M_sol,D_sol,G_sol,F_sol,C_sol,S_sol,W_sol,phi_sol,H_sol,A_sol] = deal(zeros(tot_space_steps+1, tot_time_steps+1));
%Replace the first column of the matrix for each variable with the initial
%condition
N_sol(:,1) = N0;
M_sol(:,1) = M0;
D_sol(:,1) = D0;
G_sol(:,1) = G0;
F_sol(:,1) = F0;
C_sol(:,1) = C0;
S_sol(:,1) = S0;
W_sol(:,1) = W0;
phi_sol(:,1) = phi0;
H_sol(:,1) = H0;
A_sol(:,1) = A0;
%Replace the last row of the matrix for each variable with the boundary
%condition
N_sol(end,:) = N_bound;
M_sol(end,:) = M_bound;
D_sol(end,:) = D_bound;
G_sol(end,:) = G_bound;
F_sol(end,:) = F_bound;
C_sol(end,:) = C_bound;
S_sol(end,:) = S_bound;
W_sol(end,:) = W_bound;
phi_sol(end,:) = phi_bound;
H_sol(end,:) = H_bound;
A_sol(end,:) = A_bound;
for j = 1:tot_time_steps
    %Determines whether BG fibres have completely dissolved within wound
         BG_fibres = phi_sol(:,j) > zeros(tot_space_steps+1,1);
           if sum(BG_fibres) == 0
               %M_sol(:,j) = 0.2*M_sol(:,j-1);
               % UNCOMMENT BELOW FOR REAPPLICATION OF BG FIBRES
                %for i = 1:tot_space_steps
                %    wound_healed = G_sol(i,j)+D_sol(i,j) >= 0.9;
                %    if wound_healed == 0
                %        phi_sol(i,j) = 1;
                %    else
                %        phi_sol(i,j) = 0;
                %        break
                %    end
                %end
           end
    %Nuemann boundary condition at x=0
    for i = 1:tot_space_steps
        N = N_sol(i,j);
        M = M_sol(i,j);
        D = D_sol(i,j);
        G = G_sol(i,j);
        F = F_sol(i,j);
        C = C_sol(i,j);
        S = S_sol(i,j);
        W = W_sol(i,j);
        phi = phi_sol(i,j);
        H = H_sol(i,j);
        A = A_sol(i,j);
        if i==1
            N_down = N_sol(2,j);
            N_up = N_sol(2,j);
            D_up = D_sol(2,j);
            D_down = D_sol(2,j);
            G_up = G_sol(2,j);
            G_down = G_sol(2,j);
            F_down = F_sol(2,j);
            F_up = F_sol(2,j);
            C_down = C_sol(2,j);
            C_up = C_sol(2,j);
            S_down = S_sol(2,j);
            S_up = S_sol(2,j);
            phi_down = phi_sol(2,j);
            phi_up = phi_sol(2,j);
            H_down = H_sol(2,j);
            H_up = H_sol(2,j);
            A_down = A_sol(2,j);
            A_up = A_sol(2,j);
        else
            N_down = N_sol(i-1,j);
            N_up = N_sol(i+1,j);
            D_up = D_sol(i+1,j);
            D_down = D_sol(i-1,j);
            G_up = G_sol(i+1,j);
            G_down = G_sol(i-1,j);
            F_down = F_sol(i-1,j);
            F_up = F_sol(i+1,j);
            C_down = C_sol(i-1,j);
            C_up = C_sol(i+1,j);
            S_down = S_sol(i-1,j);
            S_up = S_sol(i+1,j);
            phi_down = phi_sol(i-1,j);
            phi_up = phi_sol(i+1,j);
            H_down = H_sol(i-1,j);
            H_up = H_sol(i+1,j);
            A_down = A_sol(i-1,j);
            A_up = A_sol(i+1,j);
        end
        %differential equations
        N_sol(i,j+1) = N + delta_t*(D_n*diffusion_step(N_down,N,N_up,delta_x) - ((mean([N_up,N])*chi_n*upwinding(S_up,S,delta_x)*(mean([D_up,D]) + sigma_n*mean([G_up,G]) + sigma_nphi*0.5*(tanh(mean([phi_up,phi])/n)+1))) - (mean([N,N_down])*chi_n*upwinding(S,S_down,delta_x)*(mean([D,D_down]) + sigma_n*mean([G,G_down]) + sigma_nphi*0.5*(tanh(mean([phi,phi_down])/n)+1))))/delta_x + N*((C+beta_nh*H)/(1+S/gamma_ncs)*(1-N/kappa_n)-delta_nc/(1+C/gamma_nc)) - delta_nw*N*W);
        F_sol(i,j+1) = F + delta_t*(D_f*diffusion_step(F_down,F,F_up,delta_x) - ((mean([F_up,F])*chi_f*upwinding(S_up,S,delta_x)*(mean([D_up,D]) + sigma_f*mean([G_up,G]) + sigma_fphi*0.5*(tanh(mean([phi_up,phi])/n)+1))) - (mean([F,F_down])*chi_f*upwinding(S,S_down,delta_x)*(mean([D,D_down]) + sigma_f*mean([G,G_down]) + sigma_fphi*0.5*(tanh(mean([phi,phi_down])/n)+1))))/delta_x + beta_fb*vasculture(D,G,rho_b,gamma_b,eta) + F*((beta_fcs*C+beta_fh*H)*heaviside(S-S_bar) + beta_fc*C*(1-F/kappa_f) - delta_fc/(1+(C+rho_f*S)/gamma_fcs)) - alpha_fmw*F*W);
        C_sol(i,j+1) = C + delta_t*(D_c*diffusion_step(C_down,C,C_up,delta_x) + beta_cb*vasculture(D,G,rho_b,gamma_b,eta) - C*(delta_cn*N+delta_cf*F) - alpha_cmw*C*W);
        S_sol(i,j+1) = S + delta_t*(D_s*diffusion_step(S_down,S,S_up,delta_x) + (N*(beta_sn+beta_ssn*S)+F*(beta_sf+beta_ssf*S)+beta_sfh/beta_sc*F*H)*(1+beta_sc*heaviside(C_bar-C)) + beta_sm*M - S*(delta_s+delta_sn*N+delta_sf*F) + beta_sw*F*W);
        M_sol(i,j+1) = M + delta_t*(-M*(F*(delta_mf+delta_ms*S)+delta_m) + W*(beta_mn*N+beta_mg*G+beta_md*D));
        G_sol(i,j+1) = G + delta_t*(C*F*(beta_gf+beta_gs*S+beta_gfh*H)/(1+W/gamma_gw)*heaviside(S-S_bar)*heaviside(C_bar-C) - G*delta_gc/(1+C/gamma_gf) - alpha_gmw*G*W);
        D_sol(i,j+1) = D - delta_t*(alpha_dmw*W*D);
        W_sol(i,j+1) = W + delta_t*(W*(beta_wc*C+beta_wm*M) + beta_wn*heaviside(N_bar-N) - W*(delta_w+delta_wf*F) - delta_wa*W*A);
        phi_sol(i,j+1) = phi - delta_t*(delta_phi*heaviside(phi));
        H_sol(i,j+1) = H + delta_t*(D_h*diffusion_step(H_down,H,H_up,delta_x) + delta_phi*vol_frac*heaviside(phi-0.01)*alpha*H_d/H_w - H*(delta_hf*F+delta_hn*N+delta_h));
        A_sol(i,j+1) = A + delta_t*(D_a*diffusion_step(A_down,A,A_up,delta_x) + delta_phi*vol_frac*heaviside(phi-0.01)*H_d/H_w - A*(delta_aw*W+delta_a));
    end
end

end