diff --git a/MakeIntraSet.m b/MakeIntraSet.m
new file mode 100644
index 0000000..c8698f5
--- /dev/null
+++ b/MakeIntraSet.m
@@ -0,0 +1,44 @@
+% function Iset = MakeIntraSet(PatName)
+% Converted from BuildIntraSet on 9/1/17
+% Iset.FVIntra(Nperm,Nbm)
+
+function Iset = MakeIntraSet(PatName)
+
+homepath = pwd;
+sz = length(homepath);
+flag = 0; ind = 0;
+while flag == 0
+    ind = ind + 1;
+    if strcmp(homepath(sz-ind),'/')
+        flag = 1; endpth = sz-ind;
+    end
+end
+TrialName = homepath(endpth+1:sz);
+
+compressin = input('Compressed input: 0 = raw   1 = compressed');
+
+Trialfile = strcat(TrialName,'_Trial');
+run(Trialfile)
+
+if compressin == 1
+    pathwrk = strcat(homepath,'/dbOutComp',TrialName);
+else
+    pathwrk = strcat(homepath,'/dbOut',TrialName);
+end
+
+Nperm = Ndrug*(Ndrug-1)/2;
+Npat = 1;
+
+Ptmp = Pat2WellSet(PatName,drug,pathwrk);              %FUNCTION: /Userfunctions/Pat2WellSet
+
+Isetmp = IntraSet(Ptmp);                                % OBJECT: WellSet
+
+Iset.FVIntra(:,:) = Isetmp.FVIntra(1,:,:);
+Iset.FVIntrast(:,:) = Isetmp.FVIntrast(1,:,:);
+
+Iset.PatientName = PatName;
+Iset.DrugName = drug;
+Iset.Npat = Npat;
+Iset.Ndrug = Ndrug;
+Iset.Nperm = Ndrug*(Ndrug-1)/2;
+
diff --git a/coupleN0.m b/coupleN0.m
new file mode 100755
index 0000000..af8bc76
--- /dev/null
+++ b/coupleN0.m
@@ -0,0 +1,50 @@
+% function omegout = coupleN0(omega,g,gext)
+% function omegout = coupleN0(omega,g,gext)
+% Similar to coupleN, but for global coupling with variable coupling coeffs
+% includes a coupling matrix g(cloop,omloop)
+
+function omegout = coupleN0(omega,g,gext)
+
+[dum,N] = size(omega);
+mnomega = 1.0;
+
+options = odeset('AbsTol',1e-3);
+
+y0 = omega;
+tspan = [0 400];
+[t,y] = ode45(@f5,tspan,y0);
+[sy,dum] = size(y);
+
+y1 = y(40,:);
+
+tspan = [0 400];
+
+[t,y] = ode45(@f5,tspan,y1);
+
+for omloop = 1:N
+    [m(omloop),b(omloop)] = linfit(t,y(:,omloop));
+    w(omloop) = mnomega + m(omloop);
+%     disp(strcat('Freq',num2str(omloop),' = ',num2str(m(omloop))))
+end
+
+omegout = m;
+
+    function yd = f5(t,y)
+
+        for omloop = 1:N
+            temp = omega(omloop);
+            for cloop = 1:N
+                temp = temp + g(cloop,omloop)*sin(y(cloop)-y(omloop)) - gext*sin(y(omloop));
+            end
+            yp(omloop) = temp;
+        end
+
+        for omloop = 1:N
+            yd(omloop,1) = yp(omloop);
+        end
+
+    end     % end f5
+
+
+end % end vandpol
+
diff --git a/coupleNdriver.m b/coupleNdriver.m
new file mode 100755
index 0000000..cea533e
--- /dev/null
+++ b/coupleNdriver.m
@@ -0,0 +1,478 @@
+%coupleNdriver.m
+% Calls: makeER, makerSF, makeSW, coupleN, coupleN0
+%            makeglobal, make2Dlattice, makecycle, maketree, histfix
+
+clear
+close(figure(7))
+
+writeavi = 0;
+example = 4;        % 1 = global     2 = distributed    3 = square lattice    4 = 1D lattice   5 = tree  6 = ER 7 = Scale-free
+% 8 = small-world  9 = hypercube
+inexcoupling = 1;       % 1 = intercoupling    2 = external coupling
+Nfac = 50;
+
+displine('Example = ',example);
+displine('Coupling = ',inexcoupling);
+displine('Nfac = ',Nfac);
+
+if writeavi
+    moviename = strcat('Couple.avi');
+    aviobj = avifile(moviename,'fps',4);
+end
+
+ch = colormap(lines);
+%ch = newcolormap('fluoro');
+ch(32,:) = [0 0 0];
+ch(33,:) = [0 0 0];
+
+% ch = colormap(jet);
+% ch(32,:) = [0 0 0];
+% ch(33,:) = [0 0 0];
+
+% h = colormap(gray);
+% z = 1:32;
+% h = zeros(64,3);
+% h(1:32,1) = z/32;
+% h(33:64,2) = (33-z)/32;
+% h(32,:) = [0 0 0];
+% h(33,:) = [0 0 0];
+
+
+
+
+if example == 1     % Global coupling
+    N = 100;        %100 256
+    width = 0.2;
+    
+    omegatemp = width*(rand(1,N)-1);
+    meanomega = mean(omegatemp);
+    omega = omegatemp - meanomega;
+    sto = std(omega);
+    
+    node = makeglobal(N);
+    
+    [synch,l2,lmax] = eigenlap(node);
+    displine('synch = ',synch)
+    [cluscoef,eicoef,clus,ei] = clustercoef(node);
+    displine('cluscoef = ',cluscoef)
+    
+    nodecouple = node;
+    
+    for loop = 1:N
+        nodecouple(loop).element = omega(loop);
+        lnk(loop) = nodecouple(loop).numlink;
+    end
+    
+    avgdegree = mean(lnk)
+    
+    
+elseif example == 2        % distributed example
+    N = 11;
+    width = 0.31;
+    
+    omegat(1) = 0.16;
+    omegat(2) = 0.24;
+    omegat(3) = 0.28;
+    omegat(4) = 0.30;
+    omegat(5) = 0.31;
+    omegat(6) = -0.16;
+    omegat(7) = -0.24;
+    omegat(8) = -0.28;
+    omegat(9) = -0.30;
+    omegat(10) = -0.31;
+    omegat(11) = 0;
+    
+    for yloop = 1:N
+        for xloop = 1:N
+            g0(yloop,xloop) = 1;
+        end
+    end
+    
+    node = makeglobal(N);
+    nodecouple = node;
+    
+    for loop = 1:N
+        nodecouple(loop).element = omegat(loop);
+        lnk(loop) = nodecouple(loop).numlink;
+    end
+    
+    avgdegree = mean(lnk);
+    displine('avgdegree =',avgdegree)
+    sto = std(omegat);
+    
+    
+elseif example ==3      % square lattice
+    
+    Row = 16;       %10
+    Col = 16;           % 10
+    N = Row*Col;
+    width = 0.2;
+    
+    omegatemp = width*(rand(1,N)-1);
+    meanomega = mean(omegatemp);
+    omega = omegatemp - meanomega;
+    sto = std(omega);
+    
+    node = make2Dlattice(Row,Col);
+    
+    [synch,l2,lmax] = eigenlap(node);
+    displine('synch = ',synch)
+    [cluscoef,eicoef,clus,ei] = clustercoef(node);
+    displine('cluscoef = ',cluscoef)
+    
+    nodecouple = node;
+    
+    for loop = 1:N
+        nodecouple(loop).element = omega(loop);
+        lnk(loop) = nodecouple(loop).numlink;
+    end
+    
+    avgdegree = mean(lnk);
+    displine('avgdegree =',avgdegree)
+    
+    
+elseif example == 4      % linear cycle
+    N = 256;
+    width = 0.2;
+    
+    omegatemp = width*(rand(1,N)-1);
+    meanomega = mean(omegatemp);
+    omega = omegatemp - meanomega;
+    sto = std(omega);
+    
+    node = makecycle(N);
+    nodecouple = node;
+    
+    for loop = 1:N
+        nodecouple(loop).element = omega(loop);
+        lnk(loop) = nodecouple(loop).numlink;
+    end
+    
+    avgdegree = mean(lnk);
+    displine('avgdegree =',avgdegree)
+    
+elseif example == 5     % make a tree
+    degree = 2;
+    depth = 7;
+    width = 0.2;
+    
+    node = maketree(degree, depth);
+    
+    [dum,N] = size(node);
+    
+    omegatemp = width*(rand(1,N)-1);
+    meanomega = mean(omegatemp);
+    omega = omegatemp - meanomega;
+    sto = std(omega);
+    
+    nodecouple = node;
+    
+    for loop = 1:N
+        nodecouple(loop).element = omega(loop);
+        lnk(loop) = nodecouple(loop).numlink;
+    end
+    
+    avgdegree = mean(lnk);
+    displine('avgdegree =',avgdegree)
+    
+    
+elseif example == 6     % Erdos graph
+    N = 256;
+    p = 0.06;
+    width = 0.2;
+    
+    omegatemp = width*(rand(1,N)-1);
+    meanomega = mean(omegatemp);
+    omega = omegatemp - meanomega;
+    sto = std(omega);
+    
+    node = makeER(N,p);
+    [synch,l2,lmax] = eigenlap(node);
+    displine('synch = ',synch)
+    [cluscoef,eicoef,clus,ei] = clustercoef(node);
+    displine('cluscoef = ',cluscoef)
+    
+    nodecouple = node;
+    
+    for loop = 1:N
+        nodecouple(loop).element = omega(loop);
+        lnk(loop) = nodecouple(loop).numlink;
+    end
+    
+    avgdegree = mean(lnk);
+    displine('avgdegree =',avgdegree)
+    
+    
+    
+elseif example == 7     % scale-free graph
+    N = 256;
+    m = 1;
+    width = 0.2;
+    
+    omegatemp = width*(rand(1,N)-1);
+    meanomega = mean(omegatemp);
+    omega = omegatemp - meanomega;
+    sto = std(omega);
+    
+    node = makeSF(N,m);
+    [synch,l2,lmax] = eigenlap(node);
+    displine('synch = ',synch)
+    [cluscoef,eicoef,clus,ei] = clustercoef(node);
+    displine('cluscoef = ',cluscoef)
+    
+    nodecouple = node;
+    
+    for loop = 1:N
+        nodecouple(loop).element = omega(loop);
+        lnk(loop) = nodecouple(loop).numlink;
+    end
+    
+    avgdegree = mean(lnk);
+    displine('avgdegree =',avgdegree)
+    
+    
+elseif example == 8     % small-world graph
+    N = 256;
+    m = 3;
+    p = 0.75;
+    width = 0.2;
+    
+    omegatemp = width*(rand(1,N)-1);
+    meanomega = mean(omegatemp);
+    omega = omegatemp - meanomega;
+    sto = std(omega);
+    
+    node = makeSW(N,m,p);
+    [synch,l2,lmax] = eigenlap(node);
+    displine('synch = ',synch)
+    [cluscoef,eicoef,clus,ei] = clustercoef(node);
+    displine('cluscoef = ',cluscoef)
+    
+    nodecouple = node;
+    
+    for loop = 1:N
+        nodecouple(loop).element = omega(loop);
+        lnk(loop) = nodecouple(loop).numlink;
+    end
+    
+    avgdegree = mean(lnk);
+    displine('avgdegree =',avgdegree)
+    
+elseif example == 9     % hypercube
+    B = 8;
+    width = 0.2;
+    N = 2^B;
+    
+    omegatemp = width*(rand(1,N)-1);
+    meanomega = mean(omegatemp);
+    omega = omegatemp - meanomega;
+    sto = std(omega);
+    
+    node = makehypercube(B);
+    [synch,l2,lmax] = eigenlap(node);
+    displine('synch = ',synch)
+    [cluscoef,eicoef,clus,ei] = clustercoef(node);
+    displine('cluscoef = ',cluscoef)
+    
+    nodecouple = node;
+    
+    for loop = 1:N
+        nodecouple(loop).element = omega(loop);
+        lnk(loop) = nodecouple(loop).numlink;
+    end
+    
+    avgdegree = mean(lnk);
+    displine('avgdegree =',avgdegree)
+    
+    
+    
+end         % end if example
+
+mnomega = 1;
+
+if writeavi == 1
+    figure(2)
+    clf
+    colormap(h)
+    caxis([-0.1 0.1])
+    colorbar
+end
+
+for facloop = 1:Nfac
+    facloop
+    tic
+    
+    if example == 1     % global
+        facoef = 0.2;
+        
+    elseif example == 2  % distributed example
+        facoef = 0.3;
+        
+    elseif example == 3     % Square lattice
+        facoef = 0.35;
+        
+    elseif example == 4     % Linear cycle
+        facoef = 1.3;
+        
+    elseif example == 5     % Tree
+        facoef = 1;
+        
+    elseif example == 6     % ER graph
+        facoef = 0.4;
+        
+    elseif example == 7     % scale-free graph
+        facoef = 0.75;
+        
+    elseif example ==8     % small-world graph
+        facoef = 0.3;
+        
+    elseif example ==9     % hypercube
+        facoef = 0.3;
+        
+    end
+    
+    fac = facoef*(16*facloop/(Nfac))*(1/avgdegree)*sto/mnomega;
+    [dum,nodesz] = size(nodecouple);
+    for nodeloop = 1:nodesz
+        [dum,linksz] = size(nodecouple(nodeloop).link);
+        for linkloop = 1:linksz
+            nodecouple(nodeloop).coupling(linkloop) = fac;
+        end
+    end
+    
+    if facloop == 1     % Print the g-matrix
+        for omloop = 1:N
+            linksz = node(omloop).numlink;
+            for cloop = 1:linksz
+                ix = nodecouple(omloop).link(cloop);
+                gtemp(omloop,ix) = 1;
+            end
+        end
+        
+        figure(100)
+        imagesc(gtemp)
+        colormap(jet)
+        yyy = 1;
+    end
+    
+    facval(facloop) = fac*avgdegree;
+    
+    if inexcoupling == 1        % internal coupling
+        
+        [omegout,yout] = coupleN(nodecouple,0);                           % Here is the subfunction call for the flow
+        
+        %keyboard
+        
+    elseif inexcoupling == 2        % external coupling
+        g = 1;
+        fac = 0.05*1.16^(30*facloop/Nfac);
+        gext = fac*(1/N)*sto/mnomega;
+        
+        omegout = coupleN0(omega,zeros(N,N),gext);       % Here is the subfunction call for the flow
+        
+    end     % end if coupling
+    
+    if example == 30         % square lattice
+        ind = 0;
+        for yloop = 1:Row
+            for xloop = 1:Col
+                ind = ind+1;
+                A(yloop,xloop) = omegout(ind);
+            end
+        end
+        
+        figure(2)
+        imagesc(A)
+        colormap(h)
+        caxis([-width/2 width/2])
+        if writeavi
+            frame = getframe(gcf);
+            aviobj = addframe(aviobj,frame);
+        end
+    end
+    
+    stdev(facloop) = std(omegout);
+    
+    [y,x] = histfix(omegout,11,-0.1,0.1);
+    %    figure(3)
+    %    plot(x,y)
+    mx(facloop) = max(y);
+    
+    xx(N*(facloop-1)+1:facloop*N) = ones(1,N)*facval(facloop);
+    yy(N*(facloop-1)+1:facloop*N) = omegout;
+    
+    for omloop = 1:N
+        yyy(facloop,omloop) = omegout(omloop);
+    end
+    xxx(facloop) = facval(facloop);
+    
+    tictoc(facloop) = toc;
+    
+    S = whos;
+    [Ssz,dum] = size(S);
+    stemp = 0;
+    for sloop = 1:Ssz
+        stemp = stemp + S(sloop).bytes;
+    end
+    bytesize(facloop) = stemp;
+    
+end     % end facloop
+
+duration = sum(tictoc)
+
+figure(4)
+plot(bytesize)
+title('Byte Size')
+
+figure(5)
+plot(tictoc)
+title('Durations')
+
+figure(6)
+plot(xx,yy,'.')
+
+fig = figure(7);
+axis([0 facoef -1.1*width/2 1.1*width/2])
+%axes('LineWidth',1.1,'FontSize',14)
+
+[syy,~] = size(yyy);
+for tloop = 1:syy-3
+    yra(tloop,:) = 0.25*(yyy(tloop,:)+yyy(tloop+1,:)+yyy(tloop+2,:)+yyy(tloop+3,:));
+    %yra(tloop,:) = 0.5*(yyy(tloop,:)+yyy(tloop+1,:));
+end
+xra = xxx(1:syy-3);
+
+
+hold on
+for omloop = 1:N
+    %plot(xxx,yyy(:,omloop),'-','Color',ch(round(34-omloop*32/N),:),'LineWidth',1.5)
+    plot(xra,yra(:,omloop),'-','Color',ch(round(64-omloop*63/N),:),'LineWidth',1.5)
+end
+hold off
+fig.Color = [1 1 1];
+
+
+xp = 1:Nfac;
+
+mins = min(stdev);
+maxs = max(stdev);
+%Responst = (1 - (stdev-mins)/(maxs-mins)).^2;
+Responst = (1 - stdev/maxs);
+
+minm = min(mx);
+maxm = max(mx);
+%Responmx = (mx-minm)/(maxm-minm);
+Responmx = (mx-minm)/(1-minm);
+figure(8)
+%plot(1./facval,1-Responst,'r',1./facval,1-Responmx,'b')
+plot(1./facval,1-Responst,'b')
+title('Respon')
+legend('std','max')
+
+Printfile3('cNdout',facval,Responst,Responmx)
+
+
+if writeavi
+    aviobj = close(aviobj);
+end
+
diff --git a/histfix.m b/histfix.m
new file mode 100755
index 0000000..64b1061
--- /dev/null
+++ b/histfix.m
@@ -0,0 +1 @@
+% function [y,x] = histfix(vec,numbins,minval,maxval)
% y is probability.
% x is bin.
% Also see histfixplot.m


function [y,x] = histfix(vec,numbins,minval,maxval)

del = (maxval - minval)/numbins;
i = 1:numbins;
x = minval + i*del - del/2;

bin = zeros(1,numbins);
maxloop = length(vec);

if (maxval - minval)==0
    y = zeros(1,numbins);
else

    for loop = 1:maxloop

        index = ceil(numbins*(vec(loop)-minval)/(maxval-minval));

        if index <= 0
            bin(1) = bin(1) + 1/maxloop;
        elseif index >= numbins
            bin(numbins) = bin(numbins) + 1/maxloop;
        else
            bin(index) = bin(index) + 1/maxloop;
        end


    end

    y = bin;

end

\ No newline at end of file
diff --git a/linfit.m b/linfit.m
new file mode 100644
index 0000000..06a7cbe
--- /dev/null
+++ b/linfit.m
@@ -0,0 +1,26 @@
+% function [m,b] = linfit(x,y)
+% see also planefit.m and multifit.m
+
+function [m,b] = linfit(xx,yy,prnt,varargin)
+
+x = rowvec(xx);
+y = rowvec(yy);
+
+
+meanx = mean(x);
+meany = mean(y);
+meanxy = mean(x.*y);
+meanx2 = mean(x.^2);
+
+m = (meanxy - meanx*meany)/(meanx2 - meanx^2);
+b = meany - m*meanx;
+
+f = m*x + b;
+
+if nargin == 3
+    if prnt == 1
+        figure(103)
+        plot(x,y,'*',x,f,'r')
+        pause(0.1)
+    end
+end
diff --git a/make2Dlattice.m b/make2Dlattice.m
new file mode 100644
index 0000000..a432e5f
--- /dev/null
+++ b/make2Dlattice.m
@@ -0,0 +1,53 @@
+% function node = make2Dlattice(N,M)
+% Creates 2D lattice graph with N rows and M columns
+
+function node = make2Dlattice(N,M)
+
+node(1).element = 1;
+node(1).numlink = 0;
+node(1).link = [];
+
+% Make nodes and internal links
+for rowloop = 1:N
+    for colloop = 2:M-1
+        index = (rowloop-1)*M + colloop;
+        
+        node = addnode(index,node);
+        node(index).element = index;
+        node(index).numlink = 4;
+        node(index).link = [index-1 index+1];
+        
+    end
+    
+    node((rowloop-1)*M+1).element = (rowloop-1)*M+1;
+    node((rowloop-1)*M+1).numlink = 4;
+    node((rowloop-1)*M+1).link = [M*rowloop (rowloop-1)*M+2];
+    
+    node(rowloop*M).element = rowloop*M;
+    node(rowloop*M).numlink = 4;
+    node(rowloop*M).link = [rowloop*M-1 (rowloop-1)*M+1];
+    
+end
+
+% Make edge  links
+for rowloop = 2:N-1
+    for colloop = 1:M
+        index = (rowloop-1)*M + colloop;
+        
+        node(index).link(3) = (rowloop-2)*M + colloop;
+        node(index).link(4) = rowloop*M + colloop;
+        
+    end
+end
+
+for colloop = 1:M
+    node(colloop).link(3) = (N-1)*M + colloop;
+    node(colloop).link(4) = M + colloop;
+    
+    node((N-1)*M + colloop).link(3) = (N-2)*M + colloop;
+    node((N-1)*M + colloop).link(4) = colloop;
+    
+end
+
+
+
diff --git a/makeER.m b/makeER.m
new file mode 100644
index 0000000..5259db7
--- /dev/null
+++ b/makeER.m
@@ -0,0 +1,32 @@
+% function node = makeER(N,p)
+% Creates an Erdos-Renyi graph of N nodes and link probability p
+% node structure is ...
+% node(1).element = node number;
+% node(1).numlink = number_of_links;
+% node(1).link = [set of linked node numbers];
+
+
+function node = makeER(N,p)
+
+% Set node structure
+node(1).element = 1;
+node(1).numlink = 0;
+node(1).link = [];
+
+% Create ER graph
+[A,degree,Lap] = Erdos(N,p);
+
+% set links
+for rowloop = 2:N
+
+    node = addnode(rowloop,node);
+    
+    for coloop = 1:rowloop-1
+        if A(rowloop, coloop) ==1
+            node = addlink(rowloop,coloop,node);
+        end
+    
+    
+    
+    end     % end coloop
+end     % end rowloop
\ No newline at end of file
diff --git a/makeLB.m b/makeLB.m
new file mode 100644
index 0000000..437f814
--- /dev/null
+++ b/makeLB.m
@@ -0,0 +1,52 @@
+% function node = makeLB(N,m,)
+% Creates a linear graph (cycle) of N nodes and 2*m neighbor links per new node
+% node structure is ...
+% node(1).element = node number;
+% node(1).numlink = number_of_links;
+% node(1).link = [set of linked node numbers];
+
+
+function node = makeLB(N,m)
+
+A = zeros(N,N);
+
+node(1).element = 1;
+node(1).numlink = 2*m;
+node(1).link = [];
+for loop = 2:N
+    node = addnode(loop,node);
+    node(loop).numlink = 2*m;
+end
+
+for loop = 1:N
+    nlinks = 0;
+        for neighloop = 1:m
+            nlinks = nlinks + 1;
+            if (loop+neighloop) <=N
+                A(loop,loop+neighloop) = 1;
+                node(loop).link(nlinks) = loop+neighloop;
+            else
+                A(loop,loop+neighloop-N) = 1;
+                node(loop).link(nlinks) = loop + neighloop-N;
+            end
+            nlinks = nlinks+1;
+            if (loop-neighloop) >=1
+                A(loop,loop-neighloop) = 1;
+                node(loop).link(nlinks) = loop-neighloop;
+            else
+                A(loop,N + loop - neighloop) = 1;
+                node(loop).link(nlinks) = N + loop - neighloop;
+            end
+        end
+end
+
+% figure(1)
+% imagesc(A)
+% size(A)
+
+
+
+
+
+
+
diff --git a/makeLN.m b/makeLN.m
new file mode 100644
index 0000000..3909c83
--- /dev/null
+++ b/makeLN.m
@@ -0,0 +1,38 @@
+% function node = makeLN(N,m,)
+% Creates a linear graph (non-cycle) of N nodes and 2*m neighbor links per new node
+% node structure is ...
+% node(1).element = node number;
+% node(1).numlink = number_of_links;
+% node(1).link = [set of linked node numbers];
+
+
+function node = makeLN(N,m)
+
+node(1).element = 1;
+node(1).numlink = 0;
+node(1).link = [];
+for loop = 2:N
+    node = addnode(loop,node);
+end
+
+
+for loop = 1:N
+    nlinks = 0;
+    for neighloop = 1:m
+        if (loop+ neighloop) <=N
+            node = addlink(loop,loop+neighloop,node);
+        end
+%         if (loop - neighloop) >=1
+%             node = addlink(loop,loop-neighloop,node);
+%         end
+
+    end
+end
+
+
+
+
+
+
+
+
diff --git a/makeSF.m b/makeSF.m
new file mode 100644
index 0000000..28bbe3b
--- /dev/null
+++ b/makeSF.m
@@ -0,0 +1,45 @@
+% function node = makeSF(N,m)
+% Creates a Scale-Free graph of N nodes and m linksper new node
+% node structure is ...
+% node(1).element = node number;
+% node(1).numlink = number_of_links;
+% node(1).link = [set of linked node numbers];
+
+function node = makeSF(N,m)
+
+% Set node initial structure
+node = makeER(2*m,1);
+numnodes = 2*m;
+
+for addloop = 2*m+1:N   % add nodes up to N
+
+    node = addnode(addloop,node);
+
+    sumdegree = 0;
+    for dloop = 1:numnodes
+        sumdegree = sumdegree + node(dloop).numlink;
+    end     % end dloop
+    for dloop = 1:numnodes
+        nodeprob(dloop) = node(dloop).numlink/sumdegree;
+    end     % end dloop
+
+    for mloop = 1:m     % add m links preferentially
+
+        flag = 1;   
+        while flag == 1
+            temp = nonuniformrand(numnodes,nodeprob);
+            Ar = ismember(node(addloop).link,temp);
+            if sum(Ar) == 0
+                 nodetolink = temp;
+                 flag = 0;
+            end
+        end
+
+        node = addlink(addloop,nodetolink,node);
+
+    end     % end mloop
+
+numnodes = numnodes + 1;
+
+end     % end node addition loop
+
diff --git a/makeSW.m b/makeSW.m
new file mode 100644
index 0000000..d3a86b5
--- /dev/null
+++ b/makeSW.m
@@ -0,0 +1,52 @@
+% function node = makeSW(N,m, p)
+% Creates a Small-world graph of N nodes and 2*m links per new node
+% p is the re-wiring probability
+% node structure is ...
+% node(1).element = node number;
+% node(1).numlink = number_of_links;
+% node(1).link = [set of linked node numbers];
+% Calls:
+%  makeLB
+%  sublink
+%  addlink
+
+
+function node = makeSW(N,m,p)
+
+% Set node initial structure
+node = makeLB(N,m);             % linear cycle with 2*m neighbors
+
+for loop = 1:N
+    %loop
+    nlink = node(loop).numlink;
+    linknum = node(loop).link;
+
+    for linkloop = 1:nlink
+        oldtarget = linknum(linkloop);
+
+        if oldtarget > loop
+
+            test = rand;
+            if test < p
+
+                node = sublink(loop,oldtarget,node);
+
+                flag = 1;
+                while flag == 1
+                    ind = round(rand*(N-1)) + 1;
+                    if ind ~= loop
+                        flag = 0;
+                        node = addlink(loop,ind,node);
+                    end
+                end     % end while
+                
+%                 drawnet(node)
+%                 keyboard
+
+            end     % end test
+
+        end     % end if oldtarget
+
+    end     % end linkloop
+
+end     % end loop
diff --git a/makecycle.m b/makecycle.m
new file mode 100644
index 0000000..e0e3de0
--- /dev/null
+++ b/makecycle.m
@@ -0,0 +1,26 @@
+% function node = makecycle(N)
+% Creates cyclic (linear) graph with N nodes
+
+function node = makecycle(N)
+
+node(1).element = 1;
+node(1).numlink = 0;
+node(1).link = [];
+
+for loop = 2:N-1
+    node = addnode(loop,node);
+    node(loop).element = loop;
+    node(loop).numlink = 2;
+    node(loop).link = [loop-1 loop+1];
+end
+
+node(1).element = 1;
+node(1).numlink = 2;
+node(1).link = [N 2];
+
+node(N).element = N;
+node(N).numlink = 2;
+node(N).link = [N-1 1];
+
+
+end
diff --git a/makeglobal.m b/makeglobal.m
new file mode 100644
index 0000000..cb28e7c
--- /dev/null
+++ b/makeglobal.m
@@ -0,0 +1,22 @@
+% function node = makeglobal(N)
+% Creates complete graph with N nodes
+
+function node = makeglobal(N)
+
+node(1).element = 1;
+node(1).numlink = 0;
+node(1).link = [];
+
+for loop = 2:N
+
+    node = addnode(loop,node);
+
+    for linkloop = 1:loop-1
+
+        node = addlink(loop,linkloop,node);
+
+    end
+
+end
+
+end
diff --git a/makehypercube.m b/makehypercube.m
new file mode 100644
index 0000000..3f49184
--- /dev/null
+++ b/makehypercube.m
@@ -0,0 +1,34 @@
+% function node = makehypercube(B)
+% B = bits  (2^B)
+
+
+function node = makehypercube(B)
+
+N = 2^B;
+
+node(1).element = 1;
+node(1).numlink = 0;
+node(1).link = [];
+node(1).val = dec2bin(0,B);
+
+for yloop = 2:N
+    node = addnode(yloop,node);
+    node(yloop).val = dec2bin(yloop-1,B);
+end
+
+A = zeros(N,N);
+for yloop = 1:N
+    for xloop = yloop+1:N
+        
+        num1 = xloop - 1;
+        num2 = yloop - 1;
+                
+        dis = hamming(num1,num2);
+        if dis == 1
+            A(yloop,xloop) = 1;
+            node = addlink(node(yloop).element,node(xloop).element,node);
+        end
+    end
+end
+
+%node = adj2node(A);
diff --git a/makesiteperc.m b/makesiteperc.m
new file mode 100644
index 0000000..8136c5e
--- /dev/null
+++ b/makesiteperc.m
@@ -0,0 +1,19 @@
+% function node = siteperc(N,p)
+% makes 2D lattice site percolation with occupancy p
+% 
+
+function node = siteperc(N,p)
+
+node = make2Dlattice(N,N);
+
+prob = 1 - p;
+Del = randintexc(floor(prob*N^2),N^2);
+
+dnum = length(Del);
+
+for loop = 1:dnum
+    
+    node = subnode(Del(loop),node);
+    
+end
+
diff --git a/maketree.m b/maketree.m
new file mode 100644
index 0000000..8127a1f
--- /dev/null
+++ b/maketree.m
@@ -0,0 +1,58 @@
+% function node = maketree(degree,depth)
+% Creates a tree of given degree and depth
+% node structure is ...
+% node(1).element = node number;
+% node(1).numlink = number_of_links;
+% node(1).link = [set of linked node numbers];
+
+
+function node = maketree(degree,depth)
+
+% First calculate size of tree
+temp = 0;
+for loop = 0:depth
+    temp = temp + degree^loop;
+end
+treesize = temp;
+disp(strcat('treesize = ',num2str(treesize)))
+
+% Set node structure
+node(1).element = 1;
+node(1).numlink = 0;
+node(1).link = [];
+
+level(1) = 0;
+
+activestack.data(1) = 1;
+activestack.top = 1;
+
+availnum = 2;   % Next available node number
+
+test = 1;
+while test >= 1
+
+    [activestack,activenode] = pop(activestack);    % take top stack element
+    if level(activenode) < depth
+
+        for loop = 1:degree
+
+            newnodenum = availnum;
+
+            activestack = put(newnodenum,activestack);      % add an element to the stack for later expansion
+            node = addnode(newnodenum,node);
+            level(newnodenum) = level(activenode) + 1;
+            node = addlink(newnodenum,activenode,node);
+
+            availnum = availnum + 1;
+
+        end     % end degreeloop
+
+    end     % end if level < depth
+
+    test = activestack.top;
+
+end     % end whileloop
+
+
+end     % end function maketree
+