T = T_pde / dt;
x = -L:dx:L;
Nx = 2 * L / dx;

D = Du / K;

Nx_u = 2 * Ru / dx; % punti griglia non infette

% Inizializzazione popolazioni
u1 = zeros(Nx + 1, T_pde + 1); % cellule non infette
u3 = zeros(Nx + 1, T_pde + 1); % cellule T
u4 = zeros(Nx + 1, T_pde + 1); % sostanza chimica

u = [zeros((Nx - Nx_u)/2, 1);
     h_uninf * ones(Nx_u + 1, 1);
     zeros((Nx - Nx_u)/2, 1)];

z = zeros(Nx + 1, 1);  % cellule T
c = zeros(Nx + 1, 1);  % concentrazione chimica

u1(:, 1) = u;
u3(:, 1) = z;
u4(:, 1) = c;

for k = 1:T
    rho = u;

    % DIFFUSIONE cellule non infette
    vicino_destra = [u(2:Nx+1); 0];
    vicino_sinistra = [0; u(1:Nx)];
    u = u + Du * (vicino_sinistra + vicino_destra - 2*u) * dt / dx^2;
    u(1) = u(2);
    u(Nx+1) = u(Nx);

    % DIFFUSIONE sostanza chimica
    vicino_destra = [c(2:Nx+1); 0];
    vicino_sinistra = [0; c(1:Nx)];
    c = c + Dc * (vicino_sinistra + vicino_destra - 2*c) * dt / dx^2;
    c(1) = c(2);
    c(Nx+1) = c(Nx);

    % CHEMOTASSI cellule T
    c_destra = [c(2:Nx+1); 0];
    z_destra = [z(2:Nx+1); 0];
    w = chiT * (c_destra - c) / dx;
    Fz = max(w, 0) .* z - max(-w, 0) .* z_destra - DT * (z_destra - z) / dx;
    z = z + (dt / dx) * ([0; Fz(1:Nx)] - [Fz(1:Nx); 0]);

    % REAZIONI LOCALI
    u = u + dt * (p * (1 - rho / K) .* u - zU * u .* z);
    z = z + dt * (-qT * z + alphaz * sum(u) * dx); % afflusso in funzione del danno

    % sostanza chimica prodotta da u
    c = c + dt * ((gammac * u) .* (cstar - c) - qc * c);

    % Salvataggio dati ogni 5 ore
    if mod(k, 1/dt) == 0
        u1(:, k*dt + 1) = u;
        u3(:, k*dt + 1) = z;
        u4(:, k*dt + 1) = c;

        figure(1)
        clf
        sgtitle('t=' + string(k*dt) + ' h')
        subplot(1, 2, 1)
        plot(x, u, 'LineWidth', 1)
        hold on
        plot(x, c * 200, 'LineWidth', 1)
        legend('u (non infette)', 'c (x200)')

        subplot(1, 2, 2)
        plot(x, z, 'LineWidth', 1)
        hold on
        plot(x, c * 10, 'LineWidth', 1)
        legend('T cell', 'c (x10)')
    end
end

t = 0:T_pde;

u1 = u1';
u3 = u3';
u4 = u4';

% Modifiche principali:
%(1)rimosso completamente i, u2, Nx_i, e la diffusione/termine reattivo associato
%(2)eliminato beta, q, zI, alphac*i, ecc
%(3)la produzione della sostanza chimica è ora solo da cellule non infette (u)
%(4)il salvataggio è ora in u1, u3, u4 (coerente con i plot)