% Griglia per ritratti di fase
[u, z] = meshgrid(linspace(0, 1.5, 20), linspace(0, 1.5, 20));
ICs = [0.1 0.1; 0.9 0.1; 0.1 1.0; 0.9 1.0];
tspan = [0 50];
Sz_vals = [0.5, 1.5];  % sotto e sopra soglia

figure;

for k = 1:2
    Sz = Sz_vals(k);
    du = p*u.*(1 - u/K) - (zeta/K)*u.*z;
    dz = -qz*z + Sz;
    mag = sqrt(du.^2 + dz.^2);

    subplot(1,2,k);
    quiver(u, z, du./mag, dz./mag, 0.5, 'k'); hold on;

    f = @(t, y) [p*y(1)*(1 - y(1)/K) - (zeta/K)*y(1)*y(2); -qz*y(2) + Sz];

    for j = 1:size(ICs,1)
        [t, Y] = ode45(f, tspan, ICs(j,:));
        plot(Y(:,1), Y(:,2), 'LineWidth', 1.5);
    end

    xlabel('u (tumore)');
    ylabel('z (cellule immunitarie)');
    title(['Ritratto di fase per S_z = ' num2str(Sz)]);
    grid on;
end
