A=1;   % Sezione del serbatoio (m^2)
a=.25; % Sezione uscita (m^2)
g=9.8; % Accelerazione di gravita' (m/s^2)


h0=0;  % Altezza iniziale (m)

T=10;  % Tempo di simulazione (s) 

dt=0.1; % Step di simulazione (s)

N=floor(T/dt); % Numero di step

u=2*ones(N,1); % Vettore flusso di ingresso (m^3/s)

h=zeros(N+1,1); % Inizializza vettore delle altezze
h(1)=h0;

q=zeros(N,1); % Inizializza vettore dei flussi di uscita

K1=a*sqrt(2*g); % Una costante usata frequentemente

for t=2:N+1,
   radh=sqrt(h(t-1));
   h(t)=h(t-1)+dt*(-K1/A*radh+u(t-1)/A);
   q(t-1)=K1*radh;
end

subplot(121)
plot(0:dt:T,h)
grid
axis([0 T 0 2])
title('Altezza h(t) (m)')
xlabel('Tempo (s)')

subplot(122)
plot(0:dt:T-dt,q,0:dt:T-dt,u)
grid
axis([0 T 0 2])
title('Flussi di ingresso e uscita (m^3/s)')
xlabel('Tempo (s)')