%tic
%clear x t q;
clear;

t = 300;
t1=t;
L=1;
N=51;

x(1:N*(N+1)/2)=0;
for d=1:N,
h(d)=1;%0.5+0.5*rand;
%vv(d)=d;
end
U=5;
for n=1:N, ee(n)=2.5*sin(1+2*pi*(n-(N+1)/2)*(sqrt(5)-1)/2); end
E=ee;


% m=9; n=11;
% f=m+n*(n-1)/2;
% x(f)=-1;

n=(N+1)/2; m=n-1;
f=m+n*(n-1)/2;

x(f)=1;

options = odeset('RelTol',1e-8,'AbsTol',1e-8,'OutputFcn','odeplot','OutputSel',[1]);
[t x] = ode45('hubbardDisord1', [0:L:t], x, options,N,U,E,h);

for H=1:t1/L,
    for m=1:N, q(H,m)=0;
        for n=m:N, f=m+n*(n-1)/2;
            q(H,m)=q(H,m)+real(x(H,f)).^2+imag(x(H,f)).^2;
        end
    end
end


    for n=1:N, 
        for m=1:n, f=m+n*(n-1)/2; R(m,n)=0;
            for H=1:t1/L,
            R(m,n)=R(m,n)+real(x(H,f)).^2+imag(x(H,f)).^2;
q(H,n)=q(H,n)+real(x(H,f)).^2+imag(x(H,f)).^2;
        end
    end
end


for H=1:t1/L, M1(H)=0; M2(H)=0; qq(H)=0; HH(H)=H*L;
for I=1:N, II(I)=I;
    M1(H)=M1(H)+q(H,I)*I;
M2(H)=M2(H)+q(H,I)*I^2;
qq(H)=qq(H)+q(H,I);
end
end
M2=M2./qq; M1=M1./qq;
DM=M2-M1.^2;

fid=fopen('dima1.dat','w');
fprintf(fid,'%8.4f',q(5,:)); 
fclose(fid);