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

Nx_u=2*Ru/dx;
Nx_i=2*Ri/dx;

u1=zeros(Nx+1,T_pde+1);
u2=u1; u3=u1; u4=u1;

u=[zeros((Nx-Nx_u)/2,1)
      h_uninf*ones(Nx_u+1,1) 
      zeros((Nx-Nx_u)/2,1)];
i=[zeros((Nx-Nx_i)/2,1)
      h_inf*ones(Nx_i+1,1) 
      zeros((Nx-Nx_i)/2,1)];
z=zeros(Nx+1,1);
c=zeros(Nx+1,1);

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

for k=1:T 
    rho = (u+i);

    % --- DIFFUSIONE TUMORALE: schema pressure-driven upwind ---
    u = update_tumor_pde_pressure(u, Du, K, p, dt, dx);
    i = update_tumor_pde_pressure(i, Di, K, 0, dt, dx);  % no crescita per infette

    % --- DIFFUSIONE CHEMOATTRATTANTE ---
    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);

    % --- MOVIMENTO IMMUNE: chemotassi upwind + diffusione ---
    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 ---
    u=u+dt*( -beta*u.*i - zU*u.*z );
    i=i+dt*( beta*u.*i - q*i - zI*i.*z );
    z=z+dt*(-qT*z+alphaz*sum(i)*dx);
    c=c+dt*((alphac*i+gammac*u).*(cstar-c)-qc*c);

    if (mod(k,1/dt)==0)
        u1(:,k*dt+1)=u;
        u2(:,k*dt+1)=i;
        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,i,'LineWidth',1)
        hold on
        plot(x,c*200,'LineWidth',1)
        legend('u','i')
        
        subplot(1,2,2)
        plot(x,z,'LineWidth',1)
        hold on
        plot(x,c*10,'LineWidth',1)
        legend('z','c')
    end
end

t=0:T_pde;
u1=u1'; u2=u2'; u3=u3'; u4=u4';
