theta_1=2*Du*tau/delta^2;
theta_2=2*Di*tau/delta^2;

theta_T=2*DT*tau/delta^2;
nu=2*chiT*tau*cstar/delta^2;
Dc_ag=Dc*tau/delta^2;

T=T_pde/tau;
griglia=L/delta; %numeri pari
griglia_uninf=cell_uninf/delta; %numero punti di griglia con cellule non infette, deve essere PARI
griglia_inf=cell_inf/delta; %numero punti di griglia con cellule infette, deve essere PARI
n_uninf=h_uninf*delta; %numeri cellule non infette in ogni nodo, deve essere INTERO
n_inf=h_inf*delta; %numeri cellule infette in ogni nodo, deve essere INTERO
%NB: al momento considero +1 per coerenza con la PDE

U=zeros(2*griglia+1,T_pde+1);
I=zeros(2*griglia+1,T_pde+1);
Tcell=zeros(2*griglia+1,T_pde+1);
c=zeros(2*griglia+1,T_pde+1);

vessel=ones(2*griglia+1,1);

% Maschera inflow: solo 40% destro del dominio
mask = zeros(2 * griglia + 1, 1);
%mask(1:round(0.4 * length(mask))) = 1;
mask(end - round(0.4 * length(mask)) + 1:end) = 1;



U(:,1)=[zeros(griglia-griglia_uninf,1)
      n_uninf*ones(2*griglia_uninf+1,1) 
      zeros((griglia-griglia_uninf),1)];
I(:,1)=[zeros((griglia-griglia_inf),1)
      n_inf*ones(2*griglia_inf+1,1) 
      zeros((griglia-griglia_inf),1)];

U_temp=U(:,1);
I_temp=I(:,1);
Tc_temp=Tcell(:,1);
c_temp=c(:,1);



for i=2:T+1
   P=U_temp+I_temp;
   U_temp=passeggiata_unod_smorzata(U_temp,theta_1);
   I_temp=passeggiata_unod_smorzata(I_temp,theta_2);
   Tc_temp = passeggiata_unod_chemotaxisdiffusion_obstacle(Tc_temp, c_temp, P, theta_T, nu, cstar, w_max);

   U_temp=nascitamorte_pressione(U_temp, P, K*delta, p*tau, 0);
   I_temp=nascitamorte_pressione(I_temp, P, K*delta, 0, q*tau);
   Tc_temp=nascitamorte_pressione(Tc_temp, P, K*delta, 0, qT*tau); %or dependent on the total number
   
  
   % Inflow T solo nella regione mascherata
    Tc_temp = inflow_masked(Tc_temp, alphaz * sum(I_temp) * delta * tau, vessel, mask);

   [U_temp,I_temp] = infezione_dued(U_temp, I_temp, beta*tau/delta);

   U_temp=T_kill(U_temp, Tc_temp, min(1,zU*tau/delta));
   I_temp=T_kill(I_temp, Tc_temp, min(1,zI*tau/delta));

   c_temp=dinamicac_unod(c_temp,Dc_ag,I_temp,alphac*tau/delta,U_temp,gammac*tau/delta,cstar,qc*tau);

    if any(isnan(U_temp)) || any(isnan(I_temp)) || any(isnan(Tc_temp))
        i
    end
    
    %we only save iterations at integer times
    t=(i-1)*tau;
    if t==floor(t)
        U(:,t+1)=U_temp;
        I(:,t+1)=I_temp;
        Tcell(:,t+1)=Tc_temp;
        c(:,t+1)=c_temp;
        
        figure(1)
        clf
        subplot(1,2,1)
        plot(U_temp/delta)
        hold on
        plot(I_temp/delta)
        subplot(1,2,2)
        plot(Tc_temp/delta)
    end
end