% Parametri generali
p = 1; K = 1; qz = 1; zeta = 1;

% Sorgente immunitaria sotto e sopra la soglia critica
Sz_vals = [0.5, 1.5];
ICs = [0.1 0.1; 0.9 0.1; 0.1 1.0; 0.9 1.0];
tspan = [0 50];

figure;

for i = 1:2
    Sz = Sz_vals(i);
    f = @(t, y) [p*y(1)*(1 - y(1)/K) - (zeta/K)*y(1)*y(2);
                 -qz*y(2) + Sz];
    subplot(2,2,i); hold on;
    subplot(2,2,i+2); hold on;

    for j = 1:size(ICs,1)
        [t, Y] = ode45(f, tspan, ICs(j,:));
        subplot(2,2,i); plot(t, Y(:,1), 'DisplayName', sprintf('u0=%.1f, z0=%.1f', ICs(j,1), ICs(j,2)));
        subplot(2,2,i+2); plot(t, Y(:,2), 'DisplayName', sprintf('u0=%.1f, z0=%.1f', ICs(j,1), ICs(j,2)));
    end

    subplot(2,2,i);
    title(['Tumore u(t), S_z = ' num2str(Sz)]);
    xlabel('Tempo'); ylabel('u(t)'); legend; grid on;

    subplot(2,2,i+2);
    title(['Immunitario z(t), S_z = ' num2str(Sz)]);
    xlabel('Tempo'); ylabel('z(t)'); legend; grid on;
end
