% Loopin ulkopuolella
% ...
maxU=0;
dt00(1)=dt_cr;
% ...
% Loopin sisällä
% ...
    [maxu,in]=max(abs(u+kdt*dt_cr*udmp+(kdt*dt_cr)^2*udd)); %kdt = kriittisen aika-askeleen säätökerroin kutsussa
    if maxU==0
        maxU=maxu;
        tm=t(end);
    end
    beta=1; % Aging parameter
    maxu=maxu*exp(-beta*(t(end)-tm)/(t1-tm));    % Aging
    uddmi=abs(udd(in));
    udmi=abs(udmp(in));
    umi=abs(u(in));
    dt=1/2*abs((sqrt(4*kdt*uddmi*maxU-4*uddmi*umi+udmi^2)-udmi)/(2*uddmi));
    if dt<dt_cr/10
        dt=dt_cr/10;
    elseif dt>dt_cr*2
        dt=dt_cr;
    end
    if maxu>maxU
        maxU=maxu;
        tm=t(end);
    end
    dt00(i+1)=dt;
%...