Matlab \\ Equazioni concatenate di Allievi

Programma per il calcolo delle sovrappressioni del colpo d'ariete. Questo codice prevede l'inserimento manuale dei dati all avvio del programma. Per creare il vettore dei tempi,come indicato nel programma, è consigliabile usare la funzione "linspace(inizio, fine, lunghezza del vettore)".

%h0=carico piezometrico iniziale
%n0=valore della funzione di variazione sezione all otturatore al tempo 0
%a=celerità
%t_off=tempo della manovra di chiusura
%vettore_t=vettore contenenti i tempi a intervallo costante

function[]=allievi()
disp('definire il carico piezometrico iniziale h0(metri)')
h0=input('h0=');
disp('definire la funzione apertura otturatore n0(inizio manovra)')
n0=input('n0=');
disp('definire la celerità di propagazione')
a=input('a=');
disp('definire tempo della manovra di chiusura(secondi)')
t_off=input('t_off=');
disp('vettore dei tempi [funzione "linspace(inizio,fine,numero valori)"]')
vettore_t=input('vettore_t=');
U=2*9.8*sqrt(h0);
dim=length(vettore_t);
n_iter=dim-1;
Al=a*U/(2*9.8*h0);
periodo=vettore_t(2)-vettore_t(1);
%il ciclo serve a definire un limite per il ciclo IF nel FOR più avanti
for i=1:dim
    pass1=(t_off-vettore_t(i));
    if pass1<periodo
       n_lim=vettore_t(i);
       break
    end
end
%disp(n_lim);
lim=1-(n_lim)/t_off;
%disp(lim);
n=1-vettore_t(1)/t_off;
z_old=(-2*Al*n+sqrt(((2*Al*n)^2)+8*Al*n0+4))/(2);
z_res=zeros(n_iter,1);
z_res(1)=z_old;
n0=n;
    for flag=1:n_iter
    n=1-vettore_t(flag+1)/t_off;
    z=(-2*Al*n+sqrt(((2*Al*n)^2)-4*(z_old^2)+8+8*Al*z_old*n0))/(2);
    z_res(flag+1)=z;
    z_old=z;
    n0=n;
    if n<lim
        break
    end
    end
flag0=flag;
for flag=flag0:n_iter
%     n=0;
%     n0=0;
    z=sqrt(2-z_old^2);
    z_res(flag+1)=z;
    z_old=z;
end
z_res=h0*(z_res).^2;
z_res=real(z_res);
vettore_t=vettore_t';
tabella=[vettore_t,z_res];
disp(tabella);