diff --git a/DynamicDrawNet.m b/DynamicDrawNet.m
index 35c6fd3..b60a811 100644
--- a/DynamicDrawNet.m
+++ b/DynamicDrawNet.m
@@ -90,7 +90,7 @@
 % Time loop
 difep = 100; lastdis = 1;
 loop = 1;
-while difep > 1/N^2.5
+while difep > 1e-5*N
     loop = loop + 1;
 
     eploop = loop;
@@ -144,7 +144,8 @@
     
     %Distemp
     Dis(loop) = Distemp;
-    difep = abs(Distemp - lastdis)/N^2;
+    %difep = abs(Distemp - lastdis)/N^2;
+    difep = abs(Distemp - lastdis)/Distemp;
     lastdis = Distemp;
     
     %keyboard
diff --git a/Erdos.m b/Erdos.m
new file mode 100644
index 0000000..c8384c6
--- /dev/null
+++ b/Erdos.m
@@ -0,0 +1,40 @@
+% function [A,degree,Lap] = Erdos(N,p)
+% Generates an Erdos-Renyi random graph of N nodes with edge probability p
+% A is the adjacency matrix
+% degree is the degree of the node
+% Lap is the Lapacian matrix
+
+
+
+function [A,degree,Lap] = Erdos(N,p)
+
+e = round(p*N*(N-1)/2);
+
+A = zeros(N,N);         %Adjacency matrix
+
+loop = 0;
+while loop ~=e
+
+%     x = round(rand*(N-1))+1;
+%     y = round(rand*(N-1))+1;
+    x = ceil(N*rand);
+    y = ceil(N*rand);
+
+    flag = A(y,x);
+    if (x~=y)&(flag==0)
+        A(x,y) = 1;
+        A(y,x) = 1;
+        loop = loop +1;
+    end
+    
+
+end
+
+degree = sum(A);
+
+Lap = -A;
+for loop = 1:N
+    Lap(loop,loop) = degree(loop);
+end
+
+
diff --git a/addlink.m b/addlink.m
new file mode 100644
index 0000000..183066a
--- /dev/null
+++ b/addlink.m
@@ -0,0 +1,13 @@
+%function node = addlink(node1,node2,node)
+function node = addlink(node1,node2,node)
+
+ind = node(node1).numlink;
+node(node1).link(ind+1) = node2;
+node(node1).numlink = ind + 1;
+
+ind = node(node2).numlink;
+node(node2).link(ind+1) = node1;
+node(node2).numlink =ind+1;
+
+end
+
diff --git a/addnode.m b/addnode.m
new file mode 100644
index 0000000..5e644ce
--- /dev/null
+++ b/addnode.m
@@ -0,0 +1,13 @@
+%function node = addnode(newnodenum,node)
+
+function node = addnode(newnodenum,node)
+
+[dum,sz] = size(node);
+
+node(sz+1).element = newnodenum;
+node(sz+1).numlink = 0;
+node(sz+1).link = [];
+
+end
+
+
diff --git a/adjacency.m b/adjacency.m
new file mode 100644
index 0000000..d9dedc6
--- /dev/null
+++ b/adjacency.m
@@ -0,0 +1,34 @@
+% function [A,degree,Lap] = adjacency(node)
+% Extract adjacency matrix, node degree and Laplacian of a graph
+% node structure is ...
+% node(1).element = node number;
+% node(1).numlink = number_of_links;
+% node(1).link = [set of linked node numbers];
+
+
+function [A,degree,Lap] = adjacency(node)
+
+[dum,N] = size(node);
+A = zeros(N,N);
+Lap = zeros(N,N);
+
+for Nloop = 1:N
+    
+    [dum,L] = size(node(Nloop).link);
+    
+    for Lloop = 1:L
+        
+        A(Nloop,node(Nloop).link(Lloop)) = 1;
+        A(node(Nloop).link(Lloop), Nloop) = 1;
+        
+    end     % end Lloop
+    
+end     % end Nloop
+
+degree = sum(A);
+
+Lap = -A;
+for loop = 1:N
+    Lap(loop,loop) = degree(loop);
+end
+
diff --git a/clustercoef.m b/clustercoef.m
new file mode 100644
index 0000000..34a2cf1
--- /dev/null
+++ b/clustercoef.m
@@ -0,0 +1,35 @@
+%function  [cluscoef,eicoef,clus,ei] = clustercoef(node)
+% cluscoef = average of 2E/k(k-1)
+% eicoef = average number of shared edges
+% clus = clustering of node
+% ei = shared edges of node
+
+function [cluscoef,eicoef,clus,ei] = clustercoef(node)
+
+[Adjac,degree,Lap] = adjacency(node);
+
+[dum,N] = size(Adjac);
+
+for iloop = 1:N
+    temp = 0;
+    for rowloop = 1:N
+        for coloop = 1:N
+            temp = temp + 0.5*Adjac(iloop,rowloop)*Adjac(rowloop,coloop)*Adjac(coloop,iloop);
+        end
+    end
+
+    ei(iloop) = temp;
+    ki = node(iloop).numlink;
+
+    if ki > 1
+        clus(iloop) = 2*ei(iloop)/ki/(ki-1);
+    else
+        clus(iloop) = 0;
+    end
+end
+
+cluscoef = mean(clus);
+eicoef = mean(ei);
+
+%keyboard
+
diff --git a/clusterstats.m b/clusterstats.m
new file mode 100644
index 0000000..bd91d60
--- /dev/null
+++ b/clusterstats.m
@@ -0,0 +1,64 @@
+% function[N,e,avgdegree,maxdegree,mindegree,numclus,meanclus,Lmax,L2,LmaxL2,meandistance,diam] = clusterstats(node)
+% Generates statistics on selected  graphs
+% N = size of network
+% e = number of edges
+% avgdegree = average degree
+% maxdegree = maximum degree
+% mindegree = minimum degree
+% numclus = number of clusters
+% meanclus = cluster coef
+% Lmax = maximum eigenvalue
+% L2 = second eigenvalue
+% LmaxL2 = ratio of Lmax to L2
+% meandistance = mean of the distances
+
+function [N,e,avgdegree,maxdegree,mindegree,numclus,meanclus,Lmax,L2,LmaxL2,meandistance,diam] = clusterstats(node)
+
+[dum,N] = size(node);
+[A,degree,Lap] = adjacency(node);
+
+e = sum(degree)/2;  % number of edges
+avgdegree = mean(degree);
+maxdegree = max(degree);
+mindegree = min(degree);
+loop2 = trace(A^2);
+
+[cluscoef,eicoef,clus,ei] = clustercoef(node);
+meanclus = cluscoef;
+dis = node2distance(node);
+meandistance = mean(mean(dis))*(N^2/(N*(N-1)));
+diam = max(max(dis));
+
+
+Lam = eig(A);
+
+LamL = eig(Lap);
+Lmax = LamL(N);
+
+% Calculate number of disconnected clusters
+dflag = 0; loop = 0; cnt = 0;
+while dflag == 0
+    loop = loop + 1;
+    if loop <=N
+        if LamL(loop) < 0.1/N
+            cnt = cnt + 1;
+        else
+            dflag = 1;
+        end
+    else
+        dflag = 1;
+    end
+end
+numclus = cnt;
+
+if numclus == 1
+    L2 = LamL(2);
+    LmaxL2 = Lmax/L2;
+else
+    L2 = 0;
+    LmaxL2 = 0;
+end
+
+
+
+
diff --git a/displine.m b/displine.m
new file mode 100644
index 0000000..7bf6090
--- /dev/null
+++ b/displine.m
@@ -0,0 +1,5 @@
+% function displine(txt,val)
+
+function displine(txt,val)
+
+disp(strcat(txt,num2str(val,2)))
diff --git a/node2distance.m b/node2distance.m
new file mode 100644
index 0000000..31f5454
--- /dev/null
+++ b/node2distance.m
@@ -0,0 +1,61 @@
+% function dist = node2distance(node,maxset,varargin)
+% node is the input network
+% dist is the distance matrix
+% maxset is the maximum distance
+% disconnected cluster distance is set to maxdist+1 or maxset
+
+
+function [dist maxdis] = node2distance(node, maxset, varargin)
+
+[sy,N] = size(node);
+
+dist = zeros(N,N);
+for vertex = 1:N
+    list = node(vertex).link;
+    
+    b = zeros(N,1);
+    vec = zeros(N,1);
+    vec(vertex) = 1;
+    b = vec;
+    
+    A = adjacency(node);
+    
+    endtest = 1;
+    bold = b;
+    index = 0;
+    while (endtest ~= -1)&(index < N)
+        index = index + 1;
+        bnew = bold + A*bold;
+        
+        dif = bnew - bold;
+        for loop = 1:N
+            if (bnew(loop) >0)&(bold(loop) == 0)
+                dist(vertex,loop) = index;
+            end
+            
+        end
+        
+        endtest = where(bnew,0);
+        bold = bnew;
+        
+    end
+    
+end
+
+maxdis = max(max(dist));
+if nargin == 2
+    setmax = maxset;
+else
+    setmax = maxdis+1;
+end
+for yloop = 1:N-1
+    for xloop = yloop+1:N
+        if dist(yloop,xloop) == 0
+            %dist(yloop,xloop) = -1;
+            %dist(xloop,yloop) = -1;
+            dist(yloop,xloop) = setmax;
+            dist(xloop,yloop) = setmax;
+
+        end
+    end
+end
diff --git a/nonuniformrand.m b/nonuniformrand.m
new file mode 100644
index 0000000..efc04d4
--- /dev/null
+++ b/nonuniformrand.m
@@ -0,0 +1,28 @@
+% function nur = nonuniformrand(N,wt)
+% wt is size N vector with sum(wt) = 1
+% output is whole number between 1 and N with prob given by wt(N)
+
+function nur = nonuniformrand(N,wt)
+
+[Y,I] = sort(wt);
+
+S(1) = wt(I(1));
+for loop = 2:N
+    S(loop) = S(loop-1) + wt(I(loop));
+end
+
+R = rand;
+
+nur = 0;
+flag = 1;
+i = 1;
+while flag == 1
+    if max(S(i),R) == R
+        i = i+1;
+    else
+        nur = I(i);
+        flag = 0;
+    end
+    
+end % end while flag
+
diff --git a/randintexc.m b/randintexc.m
new file mode 100644
index 0000000..650dbf2
--- /dev/null
+++ b/randintexc.m
@@ -0,0 +1,53 @@
+% function y = randintexc(N,M,exclus,varargin)
+% returns N random integers between 1 and M exclussively
+% also excludes values in array exclus
+
+function y = randintexc(N,M,exclus,varargin)
+
+if N > M
+    disp('Error: Too many numbers')
+    y = -1;
+    return
+end
+if nargin == 3
+    if N > M - length(exclus)
+        disp('Error: Too many numbers')
+        y = -1;
+        return
+    end
+end
+
+r(1) = ceil(M*rand);
+
+for nloop = 2:N
+    
+    testflag = 0;
+    while testflag == 0
+        
+        temp = ceil(M*rand);
+        flag = 0;
+        for cntloop = 1:nloop-1
+            if temp == r(cntloop)
+                flag = 1;
+            end
+        end
+        if nargin == 3
+            for loop = 1:length(exclus)
+                if temp == exclus(loop)
+                    flag = 1;
+                end
+            end
+        end
+        
+        if flag == 0
+            r(nloop) = temp;
+            testflag = 1;
+        end
+        
+    end % end while
+    
+    
+end
+
+y = r;
+
diff --git a/removezerok.m b/removezerok.m
new file mode 100644
index 0000000..7ce7859
--- /dev/null
+++ b/removezerok.m
@@ -0,0 +1,26 @@
+
+
+function newnode = removezerok(node)
+
+tempnode = node;
+
+flag = 1;
+while flag == 1
+    
+    [dum sz] = size(tempnode);
+    
+    tempflag = 1; ind = 1;
+    while (tempflag == 1)&&(ind <=sz)
+        if tempnode(ind).numlink == 0
+            tempnode = removenode(ind,tempnode);
+            tempflag = 0;
+        end
+        ind = ind + 1;
+    end
+    if tempflag == 1
+        flag = 0;
+    end
+        
+end
+
+newnode = tempnode;
\ No newline at end of file
diff --git a/where.m b/where.m
new file mode 100644
index 0000000..ababad9
--- /dev/null
+++ b/where.m
@@ -0,0 +1,32 @@
+% function [w,N] = where(invec,num)
+% finds where num is located in 1D vector invec
+% returns all N instances
+% See also where2.m for 2D matrix
+
+function [w,N] = where(invec,num)
+eps = 1e-3;
+
+[sy,sx] = size(invec);
+if sx == 1
+    invec = invec';
+end
+[dum,sz] = size(invec);
+
+w = -1;
+ind = 0;
+for loop = 1:sz
+        tst = abs(invec(loop)-num);
+    if tst < eps
+        ind = ind+1;
+        w(ind) = loop;
+    end
+end
+N = ind;
+
+[dum,Ntemp] = size(w);
+if (Ntemp ==1&w==-1)
+    N = 0;
+end
+
+
+