diff --git a/DoublePend.m b/DoublePend.m
new file mode 100644
index 0000000..c48a05a
--- /dev/null
+++ b/DoublePend.m
@@ -0,0 +1,162 @@
+function DoublePend
+format compact
+c(1) = 'r';
+c(2) = 'b';
+c(3) = 'g';
+c(4) = 'k';
+c(5) = 'c';
+c(6) = 'm';
+c(7) = 'y';
+
+h = colormap(lines);
+
+E = 0.8;        %  1.0  1.2;    % Kinetic Energy
+
+h2 = figure(2);
+h2.Position = [11 -31 560 420];
+dum = set(h2);
+
+h3 = figure(3);
+h3.Position = [10 381 560 420];
+dum = set(h3);
+
+h5 = figure(5);
+h5.Position = [229 1 1008 804];
+dum = set(h5);
+clf
+hold on
+
+
+repnum = 128;                % 32
+mulnum = 64/repnum;
+for reploop = 1:repnum
+    
+    clear xvar yvar vy vx px py
+    
+    vx0 = 2*(0.5 - rand)*sqrt(E);
+    vy0 = -vx0 + sign(randn)*sqrt(2*E - vx0.^2);
+    
+    v = (vx0 + vy0);
+    KE = 0.5*(vx0.^2 + v.^2);
+    
+    x0 = 0.0;     % 0.1
+    y0 = -0.0;    % -0.2
+    
+    y0 = [x0 y0 vx0 vy0];
+    
+    tspan = [1 3000];            % 500
+    
+    options = odeset('RelTol',1e-9);
+    
+    [t,y] = ode45(@f1,tspan,y0,options);
+    
+    
+    siz = size(t);
+%     y1 = mod(y(1:siz,1) + pi,2*pi)-pi;  % angle-1
+%     y2 = mod(y(1:siz,2) + pi,2*pi)-pi;  % angle-2
+%     y3 = mod(y(1:siz,3) + pi,2*pi)-pi;  % angular speed-1
+%     y4 = mod(y(1:siz,4) + pi,2*pi)-pi;  % antular speed-2
+    y1 = y(1:siz,1);  % angle-1
+    y2 = y(1:siz,2);  % angle-2
+    y3 = y(1:siz,3);
+    y4 = y(1:siz,4);
+    
+    figure(1)
+    plot(t(1:siz),y1,t(1:siz),y2)
+    
+    figure(2)
+    plot(t(1:siz),y3,t(1:siz),y4)
+    xlabel('Time')
+    legend('speed')
+    
+    figure(3)
+    plot(y1,y2)
+    xlabel('x')
+    ylabel('y')
+    
+    KinEnergy = 0.5*y3.^2 + 0.5*(y3.^2 + y4.^2 + 2*y3.*y4.*cos(y2-y1));
+    PotEnergy = - 2*cos(y1) - cos(y2);
+    Energy = KinEnergy + PotEnergy;
+    figure(4)
+    plot(t(1:siz),KinEnergy,t(1:siz),PotEnergy,t(1:siz),Energy)
+    title('Energy')
+    legend('Kin','Pot','E')
+    
+    % Power Spectrum
+    Pow = 0;
+    if Pow == 1
+        yf = y(:,1);
+        F = fft(yf);
+        Pow = F.*conj(F);
+        
+        figure(4)
+        plot(Pow)
+    end
+    
+    trig = y2;  % Trigger for crossing Poincare section: try y1 and switch xvar/yvar and vx/vy
+    
+    %First Return Map
+    Fst = 1;
+    if Fst == 1
+        cnt = 0;
+        last = trig(1);
+        for loop = 2:siz
+            if (last < 0)&&(trig(loop) > 0)
+                cnt = cnt+1;
+                del1 = -trig(loop-1)./(trig(loop) - trig(loop-1));
+                
+                vx(cnt) = y3(loop-1) + del1*(y3(loop)-y3(loop-1));
+                vy(cnt) = y4(loop-1) + del1*(y4(loop)-y4(loop-1));
+                xvar(cnt) = y1(loop-1) + del1*(y1(loop)-y1(loop-1));
+                yvar(cnt) = y2(loop-1) + del1*(y2(loop)-y2(loop-1));
+                
+                px(cnt) = 2*vx(cnt) + vy(cnt)*cos(yvar(cnt)-xvar(cnt));
+                py(cnt) = vy(cnt) + vx(cnt)*cos(yvar(cnt)-xvar(cnt));
+                
+                last = trig(loop);
+            else
+                last = trig(loop);
+            end
+        end
+        
+        figure(5)
+        plot(xvar,px,'o','MarkerSize',2,'Color',h(floor(0.9*mulnum*reploop)+1,:),'MarkerFaceColor',h(floor(0.9*mulnum*reploop)+1,:))
+        xlabel('y')
+        ylabel('vy')
+    end
+    
+    pause(0.001)
+    
+end     % end reploop
+
+figure(5)
+hold off
+print -dtiff DoublePendPrint
+
+
+% Model function
+    function dy = f1(t,y)
+        
+        A = y(4).^2.*sin(y(2)-y(1));
+        B = -2*sin(y(1));
+        C = y(3).^2.*sin(y(2)-y(1)).*cos(y(2)-y(1));
+        D = sin(y(2)).*cos(y(2)-y(1));
+        EE = 2 - (cos(y(2)-y(1))).^2;
+        
+        FF = y(4).^2.*sin(y(2)-y(1)).*cos(y(2)-y(1));
+        G = -2*sin(y(1)).*cos(y(2)-y(1));
+        H = 2*y(3).^2.*sin(y(2)-y(1));
+        I = 2*sin(y(2));
+        JJ = (cos(y(2)-y(1))).^2 - 2;
+        
+        dy = zeros(4,1);
+        dy(1) = y(3);
+        dy(2) = y(4);
+        dy(3) = (A + B + C + D)./EE;
+        dy(4) = (FF + G + H + I)./JJ;
+        
+    end     % end f5
+
+
+end % end ltest
+
diff --git a/Heiles.m b/Heiles.m
new file mode 100644
index 0000000..6343183
--- /dev/null
+++ b/Heiles.m
@@ -0,0 +1,161 @@
+function Heiles
+
+c(1) = 'r';
+c(2) = 'b';
+c(3) = 'g';
+c(4) = 'k';
+c(5) = 'c';
+c(6) = 'm';
+c(7) = 'y';
+
+h = colormap(lines);
+
+
+E = 0.1;              % 1        0.1
+
+epsE = 1;               % 0.3425   1
+
+prms = sqrt(E);
+pmax = sqrt(2*E);
+
+% Potential Function
+for xloop = 1:100
+    x = -2 + 4*xloop/100;
+    for yloop = 1:100
+        y = -2 + 4*yloop/100;
+
+        V(yloop,xloop) = 0.5*x^2 + 0.5*y^2 + epsE*(x^2*y - 0.33333*y^3);
+
+    end
+end
+figure(6)
+colormap(gray)
+imagesc(V)
+hold on
+contour(V,30,'LineColor','k')
+hold off
+caxis([-0.5 1])
+colorbar
+
+%keyboard
+
+h2 = figure(2);
+h2.Position = [11 -31 560 420];
+dum = set(h2);
+
+h3 = figure(3);
+h3.Position = [10 381 560 420];
+dum = set(h3);
+
+
+h5 = figure(5);
+h5.Position = [229 1 1008 804];
+dum = set(h5);
+clf
+hold on
+repnum = 256;                % 32
+mulnum = 64/repnum;
+for reploop = 1:repnum
+    
+    clear yvar py
+    
+%     px = 2*(rand-0.499)*pmax;
+%     yp1 = 2*(rand-0.499)*sqrt(2*(E-px^2/2));
+%     py = 2*(rand-0.499)*sqrt(2*(E-px^2/2-yp1^2/2));
+%     xp1 = sign(rand-0.5)*real(sqrt(2*(E-px^2/2-yp1^2/2-py^2/2)));
+    
+    %Etest = px^2/2 + py^2/2 + xp1^2/2 + yp1^2/2
+    
+    px = 2*(rand-0.499)*pmax;
+    py = sign(rand-0.5)*real(sqrt(2*(E-px^2/2)));
+    xp1 = 0;
+    yp1 = 0;
+    
+    y0 = [xp1 yp1 px py];
+    
+    
+    
+    tspan = [1 2000];            % 500
+
+%     figure(1)
+    %options = odeset('OutputFcn',@odephas3);
+    options = odeset('RelTol',1e-9);
+
+    [t,y] = ode45(@f1,tspan,y0,options);
+
+    siz = size(t);
+    y1 = y(:,1);
+    y2 = y(:,2);
+    y3 = y(:,3);
+    y4 = y(:,4);
+
+    figure(2)
+    plot(t(1:siz),y(1:siz,1))
+    xlabel('Time')
+    legend('speed')
+
+    figure(3)
+    plot(y(1:siz,1),y(1:siz,2))
+    xlabel('x')
+    ylabel('y')
+
+
+    % Power Spectrum
+    Pow = 0;
+    if Pow == 1
+        yf = y(:,1);
+        F = fft(yf);
+        Pow = F.*conj(F);
+
+        figure(4)
+        plot(Pow)
+    end
+
+
+    %First Return Map
+    Fst = 1;
+    if Fst == 1
+        cnt = 0;
+        last = y1(1);
+        for loop = 2:siz
+            if (last < 0)&(y1(loop) > 0)
+                cnt = cnt+1;
+                del1 = -y1(loop-1)./(y1(loop) - y1(loop-1));
+                py(cnt) = y4(loop-1) + del1*(y4(loop)-y4(loop-1));
+                yvar(cnt) = y2(loop-1) + del1*(y2(loop)-y2(loop-1));
+                last = y1(loop);
+            else
+                last = y1(loop);
+            end
+        end
+
+        figure(5)
+        plot(yvar,py,'o','MarkerSize',2,'Color',h(floor(0.9*mulnum*reploop)+1,:),'MarkerFaceColor',h(floor(0.9*mulnum*reploop)+1,:))
+        xlabel('y')
+        ylabel('py')
+    end
+    
+    pause(0.001)
+
+end     % end reploop
+
+figure(5)
+hold off
+print -dtiff HeilesPrint
+
+
+% Model function
+    function dy = f1(t,y)
+
+        B = 0.00;          %0.002
+        dy = zeros(4,1);
+        dy(1) = y(3);
+        dy(2) = y(4);
+        dy(3) = -y(1) - epsE*(2*y(1)*y(2) + B*y(3));
+        dy(4) = -y(2) - epsE*(y(1)^2 - y(2)^2 + B*y(4));
+
+    end     % end f5
+
+
+end % end ltest
+
diff --git a/NaK.m b/NaK.m
new file mode 100644
index 0000000..cf5aaee
--- /dev/null
+++ b/NaK.m
@@ -0,0 +1,89 @@
+
+function NaK
+
+I = 50;     % 10   
+gL = 8;     % 8
+EL = -78;   % -80   -78
+gNa = 20;    % 20
+ENa = 60;   % 60
+gK = 10;     % 10
+EK = -90;   %-90
+C = 1;      % 1
+Vn = -25;   % -25   -45
+kn = 5;     % 5
+Vm = -20;   % -20
+km = 15;     % 15
+tau = 1;    % 1  0.152
+
+%I = 12; EL = -78; Vn = -42.5; tau = 1;                    % pg. 117  supercritical Hopf bifurcation
+%I = 4.53; EL = -80; Vn = -25; tau = 1;                  % pg. 113  Homoclinic
+I = 5; EL = -80; Vn = -25; tau = 0.152;                 % pg. 110  Bistability Try I = 4 vs. I = 5  y0=[-20,0.4],[-70,0.1]
+%I=45;gL=1;gNa=4;gK=4;Vm=-30;km=7;EL=-78;Vn=-45;tau=1;    % pg. 176  subcritical Hopf
+
+
+y0 = [-61,0.001];
+%y0 = [-20  0.4];
+%y0 = [10  0.81];
+%y0 = [-50  0.25];
+
+
+
+tspan = [0 100];
+
+figure(1)
+options = odeset('OutputFcn',@odephas2);
+[t,y] = ode45(@f5,tspan,y0,options);
+
+figure(2)
+plot(t,y(:,1),t,y(:,2))
+
+
+figure(3)
+plot(y(:,1),y(:,2))
+set(gcf, 'color', 'white')
+
+for loop = 1:100
+    xp(loop) = -80 + 100*loop/100;
+    ninf = 1./(1 + exp((Vn - xp(loop))/kn));
+    minf = 1./(1 + exp((Vm - xp(loop))/km));
+    yp(loop) = ninf;
+    yyp(loop) = (I - gL*(xp(loop) - EL) - gNa*minf*(xp(loop) - ENa))/(gK*(xp(loop) - EK));
+end
+
+hold on
+plot(xp,yp,'r',xp,yyp,'k')
+hold off
+
+printfile = 0;
+
+if printfile == 1
+    [sz,dum] = size(t);
+    ty1 = y(1:300,1)';
+    ty2 = y(1:300,2)';
+    tyV = y(1:sz,1);
+    Printfile3('Nullclines.txt',xp,yp,yyp);
+    Printfile2('Trajectory.txt',ty1,ty2);
+    Printfile2('TimeSeries.txt',t',tyV');
+end
+    
+
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+    function yd = f5(t,y)
+        
+        ninf = 1./(1 + exp((Vn - y(1))/kn));
+        minf = 1./(1 + exp((Vm - y(1))/km));
+        
+        yp(1) = I/C - gL*(y(1) - EL)/C - gNa*minf*(y(1) - ENa)/C - gK*y(2).*(y(1) - EK)/C;
+        yp(2) = (ninf - y(2))./tau;
+        
+        
+        
+        yd = [yp(1);yp(2)];
+        
+    end     % end f5
+
+end
+
diff --git a/PredPrey.m b/PredPrey.m
new file mode 100644
index 0000000..c9a6e34
--- /dev/null
+++ b/PredPrey.m
@@ -0,0 +1,47 @@
+% PredPrey.m
+% 9-3-11
+%
+
+
+function PredPrey
+
+alpha = 3;
+beta = 1;
+gamma = 3;     %-1
+delta = 1;
+epsilon = 0;
+zeta = 0;       %+1
+
+xinf = delta/epsilon;
+yinf = alpha/beta;
+eps = .1;
+
+options = odeset('RelTol',1e-6,'AbsTol',1e-7);
+
+
+%y0 = [(xinf + eps) (yinf + 0.0)];
+
+y0 = [2 2];
+
+tspan = [0 200];
+
+[t,y] = ode45(@f5,tspan,y0);
+
+figure(1)
+plot(t,y(:,1),t,y(:,2))
+
+figure(2)
+plot(y(:,1),y(:,2))
+axis([0 10 0 10])
+
+    function yd = f5(t,y)
+
+        yp(1) = y(1)*(alpha - beta*y(2) + epsilon*y(1));
+        yp(2) = -y(2) *(gamma - delta*y(1) + zeta*y(2));
+        
+        yd = [yp(1);yp(2)];
+
+    end     % end f5
+
+end
+
diff --git a/avgeigs.m b/avgeigs.m
new file mode 100644
index 0000000..d6b2fc0
--- /dev/null
+++ b/avgeigs.m
@@ -0,0 +1,50 @@
+
+N = 50;
+
+Nrep = 100;
+eigvalER = zeros(1,N);
+eigvalSF = zeros(1,N);
+eigvalSW = zeros(1,N);
+eigER = zeros(1,N);
+eigSF = zeros(1,N);
+eigSW = zeros(1,N);
+for loop = 1:Nrep
+    
+    nodeER = makeER(N,0.08);
+    nodeSF = makeSF(N,2);
+    nodeSW = makeSW(N,2,0.1);
+    
+    [AER,degree,LapER] = adjacency(nodeER);
+    [VER,DER] = eig(LapER);
+    [vER,dER] = eig(AER);
+    [ASF,degree,LapSF] = adjacency(nodeSF);
+    [VSF,DSF] = eig(LapSF);
+    [sSF,dSF] = eig(ASF);
+    [ASW,degree,LapSW] = adjacency(nodeSW);
+    [VSW,DSW] = eig(LapSW);
+    [vSW,dSW] = eig(ASW);
+    for nloop = 1:N
+        eigvalER(nloop) = eigvalER(nloop) + DER(nloop,nloop)/Nrep;
+        eigvalSF(nloop) = eigvalSF(nloop) + DSF(nloop,nloop)/Nrep;
+        eigvalSW(nloop) = eigvalSW(nloop) + DSW(nloop,nloop)/Nrep;
+        eigER(nloop) = eigER(nloop) + dER(nloop,nloop)/Nrep;
+        eigSF(nloop) = eigSF(nloop) + dSF(nloop,nloop)/Nrep;
+        eigSW(nloop) = eigSW(nloop) + dSW(nloop,nloop)/Nrep;
+    end
+    
+end
+
+x = 1:N;
+figure(1)
+plot(x,eigvalER,x,eigvalSF,x,eigvalSW)
+legend('ER','SF','SW')
+title('Graph Laplacian')
+Printfile4('GLeig',x,eigvalER,eigvalSF,eigvalSW)
+
+figure(2)
+plot(x,eigER,x,eigSF,x,eigSW)
+legend('ER','SF','SW')
+title('Adjacency')
+Printfile4('Aeig',x,eigER,eigSF,eigSW)
+
+
diff --git a/continued.m b/continued.m
new file mode 100644
index 0000000..2636efa
--- /dev/null
+++ b/continued.m
@@ -0,0 +1,20 @@
+% a = continued(r,N)
+% continued fraction of rational number r with N terms up to N = 38
+% see also convergent.m
+
+function a = continued(r,N)
+
+rtmp = r;
+
+for loop = 1:N
+    
+    tmp = fix(rtmp);
+    if abs(rtmp - tmp) > 0
+        rtmp = 1./(rtmp - tmp);
+    else
+        rtmp = 0;
+    end
+    
+    a(loop) = tmp;
+    
+end
\ No newline at end of file
diff --git a/convergent.m b/convergent.m
new file mode 100644
index 0000000..496edcd
--- /dev/null
+++ b/convergent.m
@@ -0,0 +1,40 @@
+%function [p q] = convergent(r,N)
+% Nth convergent of real number r (up to N = 38)
+% p and q are relatively prime
+%
+% see also continued.m
+
+
+function [p q] = convergent(r,Ntmp)
+
+N = Ntmp + 1;
+a = continued(r,N);
+
+num = 1;
+denom = a(N);
+
+flag = 1;
+while flag == 1
+    if denom == 0
+        N = N-1;
+        denom = a(N);
+    else
+        flag = 0;
+    end
+end
+    
+for loop = 1:N-1
+    
+    if(loop<(N-1))
+        num = a(N-loop)*denom + num;
+        tmp = denom;
+        denom = num;
+        num = tmp;
+    else
+        num = a(N-loop)*denom + num;
+    end
+end
+
+p = num;
+q = denom;
+
diff --git a/coupledrossler.m b/coupledrossler.m
new file mode 100644
index 0000000..27a2a4a
--- /dev/null
+++ b/coupledrossler.m
@@ -0,0 +1,82 @@
+% function rossler
+
+function y = coupledrossler
+
+a = 0.25;        % 0.15 (defined phase)     0.3 (chaotic phase)
+b = 0.4;         % 0.4
+c = 8.5;       % 6 8.0 
+d = 0.97;
+g = 0.1;
+
+y0 = [1 1 0 0 1 0];
+tspan = [1 30];
+[t,y] = ode45(@f5,tspan,y0);
+[sz dum] = size(y);
+
+y0 =[y(sz,1) y(sz,2) y(sz,3) y(sz,4) y(sz,5) y(sz,6)];
+tspan = [1 200];
+options = odeset('OutputFcn',@odephas3);
+[t,y] = ode45(@f5,tspan,y0,options);
+
+figure(2)
+plot(t,y(:,1),t,y(:,2),t,y(:,3))
+
+figure(3)
+plot(y(:,1),y(:,6))
+title('x1-z2')
+
+figure(4)
+plot(y(:,2),y(:,6))
+title('y1-z2')
+
+figure(5)
+plot(y(:,1),y(:,5))
+title('x1-y2')
+
+theta1 = atan2(y(:,2),y(:,1));
+theta2 = atan2(y(:,5),y(:,4));
+
+theta12 = atan2(y(:,2),y(:,4));
+theta21 = atan2(y(:,5),y(:,1));
+
+figure(6)
+plot(theta1,theta2)
+title('theta1 vs. theta2')
+
+figure(7)
+plot(theta21-theta12)
+title('theta21-theta12')
+
+phse = theta1 - theta2;
+figure(8)
+[temp bin] = maskbilevel(phse,-pi,pi,0,0);
+plot(temp)
+title('theta1-theta2')
+
+Metric1 = sqrt(sum((temp.*bin).^2))/sqrt(sum(bin.^2))
+
+
+
+d = y(:,1);
+e = y(:,2);
+f = y(:,3);
+Printfile4('lout',t',d',e',f');
+
+
+    function yd = f5(t,y)
+
+        yp(1) = -d*y(2) - y(3) + g*(y(4) - y(1));
+        yp(2) = d*y(1) + a* y(2) + g*(y(5) - y(2));
+        yp(3) = b + y(3)*(y(1) - c) + g*(y(6) - y(3));
+
+        yp(4) = -1*d*y(5) - y(6) + g*(y(1) - y(4));
+        yp(5) = 1*d*y(4) + a* y(5) + g*(y(2) - y(5));
+        yp(6) = b + y(6)*(y(4) - c) + g*(y(3) - y(6));
+
+        yd = [yp(1);yp(2);yp(3);yp(4);yp(5);yp(6)];
+
+    end     % end f5
+
+
+end % end ltest
+
diff --git a/diffusiondriver.m b/diffusiondriver.m
new file mode 100644
index 0000000..fab31c1
--- /dev/null
+++ b/diffusiondriver.m
@@ -0,0 +1,158 @@
+% diffusiondriver.m
+% 3-24-12
+
+clear
+
+close all
+
+N = 128;
+p = 0.05;
+m = 4;
+beta = 0.01;
+
+
+node = makeER(N,0.06);
+%node = makeSF(N,m);
+%node = makeSW(N,m,0.1);
+
+
+% A = adjacency(node);
+% A = ham2adj(N);
+% node = adj2node(A);
+
+[N,e,avgdegree,maxdegree,mindegree,numclus,meanclus,Lmax,L2,LmaxL2] = clusterstats(node);
+
+disp(' ')
+displine('Number of nodes = ',N)
+disp(strcat('Number of edges = ',num2str(e)))
+disp(strcat('Mean degree = ',num2str(avgdegree)))
+displine('Maximum degree = ',maxdegree)
+disp(strcat('Number of clusters = ',num2str(numclus)))
+disp(strcat('mean cluster coefficient = ',num2str(meanclus)))
+disp(' ')
+disp(strcat('Lmax = ',num2str(Lmax)))
+disp(strcat('L2 = ',num2str(L2)))
+disp(strcat('Lmax/L2 = ',num2str(LmaxL2)))
+disp(' ')
+
+ [A,degree,Lap] = adjacency(node);
+ 
+ [V,D] = eig(Lap);
+ 
+ for loop = 1:N
+     eigval(loop) = D(loop,loop);
+ end
+ 
+ figure(1)
+ plot(eigval)
+ title('Eigenvalues')
+ 
+ % initial values
+ c = zeros(N,1);
+ c(1) = 1;
+ 
+ % eigvec decomposition
+ 
+ for eigloop = 1:N
+     Vtemp = V(:,eigloop);
+     v(eigloop) = sum(c.*Vtemp);
+ end
+ 
+ % time loop
+ Ntime = 100;
+ for timeloop = 1:Ntime       % 200
+     
+     for nodeloop = 1:N
+         
+         temp = 0;
+         for eigloop = 1:N 
+             
+             temp = temp + V(nodeloop,eigloop)*v(eigloop)*exp(-eigval(eigloop)*beta*(timeloop-1));
+                  
+         end    % end eigloop
+         
+                       concentration(timeloop,nodeloop) = temp;
+
+     end    % endnodeloop
+     
+
+ end % end timeloop
+ 
+ figure(2)
+ imagesc(real(log(concentration)))
+ colormap(jet)
+ colorbar
+ caxis([-10 0])
+ title('Log Concentrations vs. time')
+ xlabel('Node Number')
+ 
+ figure(3)
+ plot(concentration(100,:))
+ title('Ending Concentrations')
+ xlabel('Node Number')
+ 
+ 
+  x = 0:Ntime-1;
+ h = colormap(jet);
+ figure(4)
+ for nodeloop = 1:N
+     rn = round(rand*63 + 1);
+     y =  concentration(:,nodeloop)+0.001;
+     semilogy(x,y,'Color',h(rn,:))
+      hold on
+ end
+ hold off
+title('Concentrations vs. time')
+ 
+ 
+ x = 0:Ntime-1;
+ h = colormap(jet);
+ figure(5)
+ for nodeloop = 1:10
+     rn = round(rand*63 + 1);
+     y =  concentration(:,nodeloop*10)+0.001;
+     %semilogy(x,y,'Color',h(rn,:),'LineWidth',1.1)
+     semilogy(x,y,'k','LineWidth',1.2)
+      hold on
+ end
+ hold off
+ set(gcf,'Color','White')
+ title('Selected Nodes: Continuous time')
+ 
+ 
+ % Now try the discrete-time-map approach
+ 
+ c0 = c;
+ dt = 1;    % 5
+ M = eye(N,N) - beta*Lap*dt;
+ 
+ for timeloop = 1:200   %40
+     
+     c = (M^timeloop)*c0;
+     
+     Con(timeloop,:) = c';
+ end
+ 
+ x = 0:1:199;   % 0:5:199
+ h = colormap(jet);
+ figure(6)
+ for nodeloop = 1:10
+     rn = round(rand*63 + 1);
+     y =  Con(:,nodeloop*10)+0.001;
+     %semilogy(x,y,'Color',h(rn,:),'LineWidth',1.1)
+     semilogy(x,y,'k','LineWidth',1.2)
+      hold on
+ end
+ hold off
+ set(gcf,'Color','White')
+ title('Selected Nodes: Discrete time')
+
+ 
+ 
+ 
+ 
+ 
+ 
+ 
+ 
+ 
\ No newline at end of file
diff --git a/flowsimple.m b/flowsimple.m
new file mode 100644
index 0000000..5e9a8fa
--- /dev/null
+++ b/flowsimple.m
@@ -0,0 +1,108 @@
+function flowsimple
+
+xrange = [-5 5];
+yrange = [-5 5];
+% xrange = [0 4];
+% yrange = [0 4];
+rngx = xrange(2) - xrange(1);
+rngy = yrange(2) - yrange(1);
+
+[X,Y] = meshgrid(xrange(1):0.1:xrange(2), yrange(1):0.1:yrange(2));
+
+smallarrow = 0;
+[x,y] = f5(X,Y);
+
+clf
+figure(1)
+for xloop = 1:10
+    xs = xrange(1) +xloop*rngx/10;
+    for yloop = 1:10
+        ys = yrange(1) +yloop*rngy/10;
+        
+        streamline(X,Y,x,y,xs,ys)
+        
+    end
+end
+hold on
+[XQ,YQ] = meshgrid(xrange(1):1:xrange(2),yrange(1):1:yrange(2));
+smallarrow = 0;
+[xq,yq] = f5(XQ,YQ);
+quiver(XQ,YQ,xq,yq,.2,'r','filled')
+hold off
+
+axis([xrange(1) xrange(2) yrange(1) yrange(2)])
+set(gcf,'Color','White')
+
+
+    function [x,y] = f5(X,Y)
+        
+%         xtmp = X*(1-X^2) - Y;
+%         ytmp = X;
+          
+%           xtmp = Y;     % Double well
+%           ytmp = -0.5*Y + X - X.^3;
+
+%           xtmp = X*(1-2*X+Y);
+%           ytmp = Y*(1+X-Y);
+
+        
+%         xtmp = X.*(1-Y.^3);        % testing dsge
+%         ytmp = Y.*(2.5-X);
+        
+%         xtmp = X.*Y + 1 - 0.5*X;
+%         ytmp = Y.*(2 - X);
+
+        xtmp = X.*(3-Y);   % Rabbit sheep
+        ytmp = Y.*(2-X);
+
+%         xtmp = Y;
+%         ytmp = -X - 0.1*exp(-abs(X)/1)*X.^2;
+
+        
+%        c = 1;       % Medio
+%        xtmp = -Y + X.*(c - Y.^2);
+%        ytmp = X + 1;
+        
+%         xtmp = X.*(1-2*X+Y);
+%         ytmp = Y.*(0.5 + 0.5*X - Y);
+        
+%         xtmp = Y+0.75*X;
+%         ytmp = X.^3-Y-X;
+        
+        
+%         xtmp =X.*(2-X-Y);
+%         ytmp  = X-Y;
+        
+        
+%         xtmp = X.*(X-Y);
+%         ytmp = Y.*(2*X-Y);
+        
+        
+%        a = 2; b = 2;      % Inward spiral
+%         xtmp = Y - a*X;
+%         ytmp = -X - b*Y;
+        
+%         xtmp = -X.^2-Y;
+%         ytmp = X;
+        
+%         xtmp = X.^2 - Y;
+%         ytmp = X - Y;
+        
+        %x = X.*(3-X-2*Y);
+        %y = Y.*(2-X-Y);
+   
+          nrm = sqrt(xtmp.^2 + ytmp.^2);
+          
+          if smallarrow == 1
+              x = xtmp./nrm;
+              y = ytmp./nrm;
+          else
+              x = xtmp;
+              y = ytmp;
+          end
+        
+    end     % end f5
+
+
+end % end flowsimple
+
diff --git a/logistic.m b/logistic.m
new file mode 100644
index 0000000..2360df4
--- /dev/null
+++ b/logistic.m
@@ -0,0 +1,40 @@
+% logistic.m
+% Logistic map
+
+clear
+
+xold = 0.231;
+maxg = 900;        % 4000 for good coverage (pause off)
+ming = 2.9;
+figure(1);
+clf
+axis([2.7 4 0 1]);
+hold on
+for gainloop = 1:maxg
+
+    g=(4.00-ming)*gainloop/maxg + ming;
+
+    for sloop = 1:100			%settle-down loop
+        xnew = g*xold*(1-xold);
+        xold = xnew;
+    end						% end settle-down loop
+
+    rep = 64;
+    ind = 0;
+    for rloop = 1:rep
+        xnew = g*xold*(1-xold);
+        xold = xnew;
+
+        ind = ind + 1;
+        x(ind) = g + (4.00-ming)*(rloop-1)/(maxg*rep);
+        y(ind) = xnew;
+
+    end		% end rloop
+
+    plot(x,y,'o','MarkerSize',2);
+     pause(0.01)  % Take pause off to watch real-time, but put maxg to 300
+
+
+end		% end gainloop
+hold off
+
diff --git a/lorenz.m b/lorenz.m
new file mode 100644
index 0000000..ff35dd0
--- /dev/null
+++ b/lorenz.m
@@ -0,0 +1,110 @@
+function lorenz
+
+a = 10;        % 10
+b = 50;         % 50
+c = 8/3;       % 8/3
+
+y0 = [20 1 1];
+%y0 = [12.8 10.955 53.455];
+%y0 = [-4.3426   -9.3985 21.389];
+
+tspan = [0 20];
+[t,yy] = ode45(@f5,tspan,y0);
+[sz,~] = size(yy);
+
+y0 = yy(sz,:);
+tspan = [0 100];
+
+options = odeset('OutputFcn',@odephas3);
+% options = odeset('OutputFcn',@odeplot);
+options = odeset('RelTol',1e-5);
+
+[t,y] = ode45(@f5,tspan,y0,options);
+
+figure(2)
+colormap(jet)
+plot(t,y(:,1),'k',t,y(:,2),'r',t,y(:,3),'b')
+legend('x','y','z')
+
+figure(20)
+colormap(jet)
+plot3(y(:,1),y(:,2),y(:,3),'k')
+set(gcf,'Color','white')
+grid on
+
+
+figure(3)
+plot(y(:,1),y(:,2))
+xlabel('x-value')
+ylabel('y-value')
+
+u = sqrt(y(:,1).^2 + y(:,2).^2);
+z = y(:,3);
+
+figure(4)
+plot(u,z)
+xlabel('u-value')
+ylabel('z-value')
+
+xp = u - 18;
+yp = z - 52;
+
+theta = atan2(yp,xp);
+signal = sin(theta);
+
+figure(5)
+plot(signal,'LineWidth',2)
+
+
+PowS = 1;
+% Power Spectrum
+if PowS == 1        % Turn this off by setting PowS = 0
+    [st dum] = size(t);
+    tshort = t(400:st,1);
+    %yshort = y(400:st,1);
+    yshort = u(400:st,1);
+    [st dum] = size(tshort);
+    dt = (max(tshort)-min(tshort))/st;
+    df = 1/st/dt;
+    fmax = 1/dt;
+    freq = df:df:fmax;
+
+    yf = yshort(:,1);
+    F = fft(yf);
+    Pow = F.*conj(F);
+
+    mid = round(st/4);
+    figure(6)
+    semilogx(freq(2:mid),Pow(2:mid))
+    title('Power Spectrum')
+    xlabel('Frequency (Hz)')
+    
+
+    [Pmax I] = max(Pow);
+    f0 = freq(I)
+    T0 = 1/f0
+%     T0nat = 2*pi/sqrt(w02)
+end
+
+
+
+
+d = y(:,1);
+e = y(:,2);
+f = y(:,3);
+Printfile4('lout',t',d',e',f')
+
+
+    function yd = f5(t,y)
+
+        yp(1) = a*(y(2) - y(1));
+        yp(2) = y(1)*(b - y(3)) - y(2);
+        yp(3) = y(1)*y(2) - c*y(3);
+
+        yd = [yp(1);yp(2);yp(3)];
+
+    end     % end f5
+
+
+end % end ltest
+
diff --git a/lotkavoltstream.m b/lotkavoltstream.m
new file mode 100644
index 0000000..f90e66a
--- /dev/null
+++ b/lotkavoltstream.m
@@ -0,0 +1,97 @@
+% function lotkavoltstream
+%
+
+
+function lotkavoltstream
+
+alpha = 1;
+beta = 1.5;
+gamma = 1.;
+delta = 1.25;
+
+figure(1)
+N = 20;
+for loop = 2:N
+    
+    xstar = gamma/delta;
+    ystar = alpha/beta;
+    
+    valx = xstar*loop/N;
+    valy = ystar*loop/N;
+    y0 = [valx valy];
+    tspan = [0 10];
+    [t,y] = ode45(@F5,tspan,y0);
+    [sz dum] = size(t);
+    y0 = [y(sz,1) y(sz,2)];
+    clear t y
+    tspan = [10 20];
+    options = odeset('RelTol',1e-7);
+    [t,y] = ode45(@F5,tspan,y0,options);
+    
+    plot(y(:,1),y(:,2),'LineWidth',1.5)
+    mxlim = 3*max(xstar,ystar);
+    axis([0 mxlim 0 mxlim])
+
+    hold on
+    
+end
+hold off
+set(gcf,'Color','white')
+
+
+figure(2)
+N = 3;
+for loop = 2:N
+    
+    xstar = gamma/delta;
+    ystar = alpha/beta;
+    
+    valx = xstar*(loop-1)/N;
+    valy = ystar*(loop-1)/N;
+    y0 = [valx valy];
+    tspan = [0 10];
+    [t,y] = ode45(@F5,tspan,y0);
+    [sz dum] = size(t);
+    y0 = [y(sz,1) y(sz,2)];
+    clear t y
+    tspan = [10 40];
+    options = odeset('RelTol',1e-7);
+    [t,y] = ode45(@F5,tspan,y0,options);
+    
+    plot(t,y(:,1),t,y(:,2),'--','LineWidth',1.5)
+    mxlim = 3*max(xstar,ystar);
+
+    hold on
+    
+end
+hold off
+set(gcf,'Color','white')
+
+
+
+
+
+%keyboard
+
+    function [x,y] = f5(X,Y)
+        
+        
+        x = X.*(alpha - beta*Y);
+        y = -Y*(gamma - delta*X);
+        
+        
+    end     % end f5
+
+    function yd = F5(t,y)
+        
+        yp(1) = y(1).*(alpha - beta*y(2));
+        yp(2) = -y(2).*(gamma - delta*y(1));
+        
+        
+        yd = [yp(1);yp(2)];
+        
+    end     % end f5
+
+
+end % end flowsimple
+
diff --git a/quasiHam.m b/quasiHam.m
new file mode 100644
index 0000000..cc5e622
--- /dev/null
+++ b/quasiHam.m
@@ -0,0 +1,145 @@
+% Quasispecies simulation.  Includes mutation with Hamming distance.
+
+function quasiHam
+
+%close all
+h = colormap(lines);
+
+randpop = 0;
+
+B = 8;
+N = 2^B;         % size of mutation space (32)
+ep = 0.1;      % average mutation rate: 0.1 to 0.01 typical
+lam = 1;     % gain of Hamming distance
+
+if randpop == 1
+    x0temp = rand(1,N);    % Initial population
+    sx = sum(x0temp);
+    x0 = x0temp/sx;
+else
+    x0 = zeros(1,N);
+    x0(1) = 0.667; x0(2) = 0.333;
+end
+
+Pop0 = sum(x0);
+
+for yloop = 1:N
+    for xloop = 1:N
+        H(yloop,xloop) = hamming(yloop,xloop);
+    end
+end
+
+%keyboard
+
+
+Qtemp = 1./(1+ep*H);    %Mutation matrix
+Qsum = sum(Qtemp,2);  
+
+% Normalize along species
+for yloop = 1:N
+    for xloop = 1:N
+        Q(yloop,xloop) = Qtemp(yloop,xloop)/Qsum(yloop);
+    end
+end
+
+xp = 1:N;
+alpha = 64;
+ftemp = exp(-lam*H(alpha,:));   % Fitness landscape
+sf = sum(ftemp);
+f = ftemp/sf;
+
+% Transition matrix
+for yloop = 1:N
+    for xloop = 1:N
+        W(yloop,xloop) = f(yloop)*Q(yloop,xloop);
+    end
+end
+
+tspan = [0 500];
+
+
+[t,x] = ode45(@quasispec,tspan,x0);
+
+Pop0
+[sz,dum] = size(t);
+Popend = sum(x(sz,:))
+
+phistar = sum(f.*x(sz,:))
+
+
+figure(1)
+plot(f,'o-')
+hold on
+figure(1)
+plot(x(sz,:),'o-r')
+hold off
+legend('fitness','population')
+
+figure(2)
+for loop = 1:N
+    semilogx(t,x(:,loop),'Color',h(1+round(loop*63/N),:))
+    hold on
+end
+hold off
+
+figure(3)
+for loop = 1:N
+    plot(t,x(:,loop),'Color',h(1 + round(loop*63/N),:),'LineWidth',1.25)
+    hold on
+end
+hold off
+set(gcf,'Color','white')
+xlabel('Time','FontSize',14)
+ylabel('Population','FontSize',14)
+hh = gca;
+set(hh,'FontSize',14)
+
+figure(4)
+for loop = 1:N
+    loglog(t,x(:,loop),'Color',h(1+round(loop*63/N),:))
+    hold on
+end
+hold off
+
+figure(5)
+for loop = 1:N
+    semilogy(t,x(:,loop),'Color',h(1+round(loop*63/N),:),'LineWidth',1.25)
+    hold on
+end
+hold off
+set(gcf,'Color','white')
+xlabel('Time','FontSize',14)
+ylabel('Population','FontSize',14)
+hh = gca;
+set(hh,'FontSize',14)
+
+%keyboard
+
+
+
+% Eigenvalues
+
+[V,D] = eig(W);
+
+max(D(:,1))
+
+figure(6)
+plot(V(:,2))
+title('eigenvalues')
+
+
+disp(' ')
+
+
+    function yd = quasispec(t,y)
+        
+        phi = sum(f'.*y);   % Average fitness of population
+
+        yd = W*y - phi*y;
+        
+
+    end     % end quasispec
+
+end     % end quasi
+
+
diff --git a/quasiSpec.m b/quasiSpec.m
new file mode 100644
index 0000000..f1a511e
--- /dev/null
+++ b/quasiSpec.m
@@ -0,0 +1,223 @@
+% Quasispecies simulation.  Includes mutation with Hamming distance.
+
+function quasiSpec
+
+%close all
+global W
+h = colormap(lines);
+
+randpop = 0;    % 0) = spike population; 1) = random population
+mutype = 0;     % 0) = Hamming; 1) = rand
+fitype = 0;     % 0) = Hamming; 1) = 2-peak; 2) = rand+gauss;  3) freq-dep
+
+B = 7;
+N = 2^B;         % size of mutation space (64)
+lam = 1;        % Hamming fitness only
+gamma = 1;    % freq-dep fitness (payoff matrix only)
+relran = 0.025;      % relative random contrib to fitness
+time_expand = 50;
+
+ep = 0.08;      % average mutation rate: 0.1 to 0.01 typical  (0.4835) (0.0290)
+
+
+%%%%% Set original population
+if randpop == 1
+    rng(0);
+    x0temp = rand(1,N);    % Initial population
+    sx = sum(x0temp);
+    x0 = x0temp/sx;
+else 
+    x0 = zeros(1,N);
+    x0(1) = 0.667; x0(2) = 0.333;
+end
+
+Pop0 = sum(x0);
+
+
+%%%%%% Set Hamming distance
+for yloop = 1:N
+    for xloop = 1:N
+        H(yloop,xloop) = hamming(yloop-1,xloop-1);
+    end
+end
+
+
+
+%%%%%%% Set Mutation matrix
+if mutype == 0
+    Qtemp = 1./(1+H/ep);    %Mutation matrix on Hamming
+    %Qtemp = exp(-H/(ep*50));
+    Qsum = sum(Qtemp,2);
+    
+    % Normalize mutation among species
+    for yloop = 1:N
+        for xloop = 1:N
+            Q(yloop,xloop) = Qtemp(yloop,xloop)/Qsum(xloop);
+        end
+    end
+    
+end
+if mutype == 1  
+    rng(0);
+    S = stochasticmatrix(N);
+    Stemp = S - diag(diag(S));
+    Qtemp = ep*Stemp;
+    sm = sum(Qtemp,2)';
+    Q = Qtemp + diag(ones(1,N) - sm);
+end
+
+
+%keyboard
+
+
+%%%%%%% Set fitness landscape
+
+if fitype == 0      % Hamming
+    x = 1:N;
+    alpha = 84;
+    ftemp = exp(-lam*H(alpha,:));   % Fitness landscape
+    sf = sum(ftemp);
+    f = ftemp/sf;
+end
+if fitype == 1      % double peak and rand
+    rng(1);
+    f = rand(1,N);
+    x = 1:N;
+    delg = 20;
+    sig1 = 1;
+    sig2 = 4;
+    g1 = gaussprob(x,(N/2 - delg),sig1);
+    g2 = 3*gaussprob(x,(N/2 + delg),sig2);
+    ftemp = relran*f + g1 + g2;
+    f = ftemp/sum(ftemp);
+end
+if fitype == 2      % rand + Gauss
+    rng(0);
+    f = rand(1,N);
+    x = 1:N;
+    ftemp = relran*f + gauss((x-N/2)/2);   % Fitness landscape
+    f = ftemp/sum(ftemp);
+end
+if fitype == 3      % frequency-dependent Hamming
+    avgdis = mean(mean(H));
+    %payoff = exp(-gamma*(H - avgdis));    % payoff matrix
+    %payoff = H.^2;
+    %payoff = ones(size(H));
+    payoff = exp(-gamma*H);
+end
+
+%keyboard
+
+% Run time evolution
+tspan = [0 1000];
+[t,x] = ode45(@quasispec,tspan,x0);
+
+Pop0
+[sz,dum] = size(t);
+Popend = sum(x(sz,:))
+
+phistar = sum(f.*x(sz,:))       % final average fitness
+
+
+figure(1)
+plot(f,'-')
+hold on
+figure(1)
+plot(x(sz,:),'r')
+hold off
+legend('fitness','population')
+
+
+figure(2)
+for loop = 1:N
+    semilogx(t,x(:,loop),'Color',h(round(loop*64/N),:),'LineWidth',1.25)
+    hold on
+end
+hold off
+set(gcf,'Color','white')
+xlabel('Time','FontSize',14)
+ylabel('Population','FontSize',14)
+hh = gca;
+set(hh,'FontSize',14)
+
+
+figure(3)
+for loop = 1:N
+    plot(t,x(:,loop),'Color',h(round(loop*64/N),:))
+    hold on
+end
+hold off
+
+figure(4)
+for loop = 1:N
+    loglog(t,x(:,loop),'Color',h(round(loop*64/N),:))
+    hold on
+end
+hold off
+
+figure(5)
+for loop = 1:N
+    semilogy(t,x(:,loop),'Color',h(round(loop*64/N),:))
+    hold on
+end
+hold off
+
+
+
+% Eigenvalues
+
+[V,D] = eig(W);
+
+max(D(:,1))
+
+figure(6)
+%semilogy(abs(V(:,1)))
+plot((V(:,1)))
+title('Eigenvector')
+
+disp(' ')
+
+
+if fitype == 1
+    xlo = N/2 - delg - 2*sig1;
+    xhi = N/2 - delg + 2*sig1;
+    fit44 = sum(f(xlo:xhi))
+    pop44 = sum(x(sz,xlo:xhi))/Popend
+    xlo = N/2 + delg - 2*sig2;
+    xhi = N/2 + delg + 2*sig2;
+    fit84 = sum(f(xlo:xhi))
+    pop84 = sum(x(sz,xlo:xhi))/Popend
+    
+end
+
+
+
+    function yd = quasispec(~,y)
+        
+        if fitype == 3      % frequency-dependent Hamming
+            for loop = 1:N
+                ftemp(loop) = sum(payoff(:,loop).*y);
+            end
+            f = time_expand*ftemp/sum(ftemp);
+        end
+        
+        
+        % Transition matrix
+        for yloop = 1:N
+            for xloop = 1:N
+                W(yloop,xloop) = f(yloop)*Q(yloop,xloop);
+            end
+        end
+        
+        phi = sum(f'.*y);   % Average fitness of population
+        
+        yd = W*y - phi*y;
+        
+        
+    end     % end quasispec
+
+
+
+
+end     % end quasiSpec
+
diff --git a/randsgn.m b/randsgn.m
new file mode 100644
index 0000000..a3155d2
--- /dev/null
+++ b/randsgn.m
@@ -0,0 +1,10 @@
+% function y = randsgn(N)
+
+
+function y = randsgn(N)
+
+%x = randbin(N,0.5);
+
+ytmp = 2*(0.5-rand(1,N));
+y = sign(ytmp);
+
diff --git a/repeq.m b/repeq.m
new file mode 100644
index 0000000..f4f30aa
--- /dev/null
+++ b/repeq.m
@@ -0,0 +1,194 @@
+% replicator equation
+% 3-10-12
+
+function repeq
+
+N = 7;
+asymm = 4;      % 1 = zero diag (replicator eqn)   2 = zero trace (autocatylitic model)  3 = winner take all   4 = hypercycle
+phi0 = 0.01;            % average fitness (positive number) damps oscillations
+fac = 1;
+
+displine('N = ',N)
+displine('asymm = ',asymm)
+displine('phi0 = ',phi0)
+
+%h = newcolormap('fluorodark');
+h = colormap(jet);
+hg = [1 1 1] - colormap(gray);
+
+tempx = rand(1,N);
+x0 = tempx/sum(tempx);% initial populations
+
+node = makeER(N,1);
+%node = makeSW(N,N/32,0.25);
+%node = makeSF(N,N/32);
+Adj = adjacency(node);
+numclus = clusternum(node);
+[N,e,n] = clusterstats(node);
+displine('number of clusters = ',numclus)
+displine('average degree = ',n)
+
+if asymm == 1
+    % Asymmetric payoff matrix with zero diagonals (replicator model)
+    A = zeros(N,N);  
+    for yloop = 1:N
+        for xloop = yloop + 1:N
+            A(yloop,xloop) = 2*(0.5 - rand);
+            A(xloop,yloop) = -A(yloop,xloop);
+        end
+    end
+    
+elseif asymm == 2
+    % Asymmetric payoff matrix with zero trace (autocatylitic model)
+    Atemp = zeros(N,N);
+    for yloop = 1:N
+        for xloop = yloop+1:N
+            Atemp(yloop,xloop) = 2*(0.5 - rand);
+            Atemp(xloop,yloop) = -Atemp(yloop,xloop);
+        end
+        Atemp(yloop,yloop) = 2*(0.5 - rand);
+    end
+    tr = trace(Atemp);
+    A = Atemp;
+    
+    for yloop = 1:N
+        A(yloop,yloop) = Atemp(yloop,yloop) - tr/N;
+    end
+    
+elseif asymm == 3
+    % Alternating positive/negative payoff matrix with zero diagonals (winner take all model)
+    A = zeros(N,N);
+    for yloop = 1:N-1
+        xloop = yloop + 1;
+        A(yloop,xloop) = odd(xloop)*abs(randn)-even(xloop)*abs(randn);
+        A(xloop,yloop) = -A(yloop,xloop);
+    end
+    
+elseif asymm == 4 
+    % hypercycle
+    A = zeros(N,N);
+    for yloop = 2:N
+        A(yloop,yloop-1) = 0.1*N + rand;
+    end
+    A(1,N) = 0.1*N + rand;
+    
+    
+end % end if asymm
+
+%keyboard
+
+E = eye(N,N);
+A = fac*(A.*Adj + E.*A);
+
+%keyboard
+
+figure(1)
+imagesc(A);
+colormap(h)
+colorbar
+title('Payoff Matrix')
+
+tspan = [1 1600];
+[t,x] = ode45(@f5,tspan,x0);
+
+[sy,sx] = size(x);
+
+% fitness
+xend = x(sy,:)';
+fit = A*xend;
+
+avgfit = fit'*xend;
+
+pop = sum(x,2);
+av = mean(mean(x));
+displine('population = ',mean(pop))
+displine('av pop =', mean(mean(x)));
+
+% How many non-zero?
+cnt = 0;
+for loop = 1:N
+    if mean(x(sy-50:sy,loop),1) > av/20
+        cnt = cnt + 1;
+    end
+end
+displine('cnt = ',cnt)
+displine('avg fitness = ',avgfit)
+xend
+fit
+A
+
+% Average max slope?
+% for loop = 1:N
+%     mx = max(x(200:sy,loop));
+%     deriv(:,loop) = diff(x(200:sy,loop));
+%     freq(loop) = mean(abs(deriv(:,loop)))/mx;
+% end
+% mnfreq = mean(freq);
+% displine('mean frequency = ',mnfreq)
+
+figure(2)
+for loop = 1:N
+    %plot(t,real(log(x(:,loop))),'Color',h(round(loop*64/N),:),'LineWidth',1.25)
+    %axis([0 1600 -10 0])
+    plot(t,x(:,loop),'Color',h(round(loop*64/N),:),'LineWidth',1.25)
+    hold on
+end
+hold off
+
+figure(3)
+for loop = 1:N
+    plot(t,real(log(x(:,loop))),'Color',h(round(loop*64/N),:),'LineWidth',1.25)
+    axis([0 1600 -10 0])
+    %plot(t,x(:,loop),'Color',h(round(loop*64/N),:),'LineWidth',1.25)
+    hold on
+end
+hold off
+
+% figure(3)
+% plot(x(:,1),x(:,N))
+% 
+% figure(4)
+% plot(t,x(:,1),t,(x(:,round(N/2))),t,x(:,N))
+
+%keyboard
+
+
+    function yd = f5(t,y)
+        
+        for iloop = 1:N
+            ftemp = 0;
+            for jloop = 1:N
+                ftemp = ftemp + A(iloop,jloop)*y(jloop);
+            end     % end jloop
+            f(iloop) = ftemp;
+        end     % end iloop
+        
+        phitemp = phi0;          % Can adjust this from 0 to 1 to stabilize (but Nth population is no longer independent)
+        for loop = 1:N
+            phitemp = phitemp + f(loop)*y(loop);
+        end
+        phi = phitemp;
+        
+        for loop = 1:N-1
+            yd(loop) = y(loop)*(f(loop) - phi);
+        end
+        
+        if abs(phi0) < 0.01             % average fitness maintained at zero
+            yd(N) = y(N)*(f(N)-phi);
+        else                                      % non-zero average fitness
+            ydtemp = 0;
+            for loop = 1:N-1
+                ydtemp = ydtemp - yd(loop);
+            end
+            yd(N) = ydtemp;
+        end
+        
+        yd = yd';
+        
+        
+    end     % end f5
+
+
+end     % end repeq
+
+
diff --git a/repmut.m b/repmut.m
new file mode 100644
index 0000000..f8d57fd
--- /dev/null
+++ b/repmut.m
@@ -0,0 +1,164 @@
+%function repmut
+% 12/08/14
+
+function repmut
+
+clear
+format compact
+
+
+f = 1;  % this is a dummy init to make f global
+N = 64;     % 64
+p = 1/sqrt(N);
+time_expand = N;
+
+mutype = 0;     % 0 = Hamming   1 = rand
+pay = 0;        % 0 = Hamming   1 = 1/sqrt(N) 
+ep = .1;      % average mutation rate: 0.1 to 0.01 typical  (0.4835)
+
+%%%%% Set original population
+x0temp = rand(1,N);    % Initial population
+sx = sum(x0temp);
+y0 = x0temp/sx;
+Pop0 = sum(y0);
+
+
+%%%%% Set Adjacency
+
+node = makeglobal(N);
+%node = makeER(N,0.1);
+%node = makeSF(N,4);
+%node = makeSW(N,4,0.5);
+[Adj,degree,Lap] = adjacency(node);
+
+
+
+%%%%%% Set Hamming distance
+for yloop = 1:N
+    for xloop = 1:N
+        H(yloop,xloop) = hamming(yloop-1,xloop-1);
+    end
+end
+
+
+
+%%%%%%% Set Mutation matrix
+if mutype == 0
+    Qtemp = 1./(1+H/ep);    %Mutation matrix on Hamming
+    %Qtemp = exp(-H/(ep*50));
+    Qsum = sum(Qtemp,2);
+    
+    % Normalize mutation among species
+    for yloop = 1:N
+        for xloop = 1:N
+            Q(yloop,xloop) = Qtemp(yloop,xloop)/Qsum(xloop);
+        end
+    end
+    
+elseif mutype == 1  
+    S = stochasticmatrix(N);
+    Stemp = S - diag(diag(S));
+    Qtemp = ep*Stemp;
+    sm = sum(Qtemp,2)';
+    Q = Qtemp + diag(ones(1,N) - sm);
+end
+
+
+figure(1)
+imagesc(Q)
+title('Mutation Matrix')
+colormap(jet)
+
+%%%%%%% Set payoff matrix
+if pay == 1
+    payoff = zeros(N,N);
+    for yloop = 1:N
+        payoff(yloop,yloop) = 1;
+        for xloop = yloop + 1:N
+            payoff(yloop,xloop) = p;
+            payoff(xloop,yloop) = p;
+            %payoff(yloop,xloop) = p*2*(0.5 - randbin(1,0.5));
+            %payoff(xloop,yloop) = payoff(yloop,xloop);
+            %payoff(xloop,yloop) = -payoff(yloop,xloop);
+        end
+    end
+elseif pay == 0
+    payoff = exp(-1*H);
+end
+
+
+figure(2)
+imagesc(payoff)
+title('Payoff Matrix')
+colormap(jet)
+
+
+% Run time evolution
+tspan = [0 1000];
+[t,x] = ode45(@quasispec,tspan,y0);
+
+Pop0
+[sz,dum] = size(t);
+Popend = sum(x(sz,:))
+
+phistar = sum(f.*x(sz,:))       % final average fitness
+
+
+
+figure(3)
+clf
+h = colormap(lines);
+for loop = 1:N
+    plot(t,x(:,loop),'Color',h(round(loop*64/N),:))
+    hold on
+end
+hold off
+
+figure(4)
+clf
+for loop = 1:N
+    semilogx(t,x(:,loop),'Color',h(round(loop*64/N),:))
+    hold on
+end
+hold off
+
+figure(5)
+clf
+for loop = 1:N
+    semilogy(t,x(:,loop),'Color',h(round(loop*64/N),:))
+    hold on
+end
+hold off
+
+%keyboard
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+    function yd = quasispec(~,y)
+        
+        for loop = 1:N
+            ftemp(loop) = sum(payoff(:,loop).*y);
+        end
+        f = time_expand*ftemp/sum(ftemp);
+        
+        
+        % Transition matrix
+        for yloop = 1:N
+            for xloop = 1:N
+                W(yloop,xloop) = f(yloop)*(Adj(yloop,xloop)*Q(yloop,xloop));
+            end
+        end
+        
+        phi = sum(f'.*y);   % Average fitness of population
+        
+        yd = W*y - phi*y;
+        
+        
+    end     % end quasispec
+
+end
+
+
diff --git a/rksimple.m b/rksimple.m
new file mode 100644
index 0000000..e3b9115
--- /dev/null
+++ b/rksimple.m
@@ -0,0 +1,53 @@
+
+function rksimple
+
+clear
+
+
+y0 = [1 0];
+%y0 = [4,4];
+tspan = [0 400];
+
+[t,y] = ode45(@f5,tspan,y0);
+
+figure(1)
+%plot(t,y(:,1),t,y(:,2),t,y(:,3))
+plot(t,y(:,1),t,y(:,2));
+legend('1','2','3')
+
+figure(2)
+plot(y(:,1),y(:,2))
+
+
+    function yd = f5(t,y)
+        
+        w = 1;
+
+        yp(1) = y(2);
+        yp(2) = -y(1) + 0.5*y(2)*(1-y(1).^2);
+        
+        %mu = 1; beta = 1; w02 = 1;
+        %yp(1) = y(2);
+        %yp(2) = mu*y(2)*(1-beta*y(1)^2)-w02*y(1);
+        
+%         yp(1) = -y(2) - y(1).^3 + 5*y(1);
+%         yp(2) = -y(2) + 4*(y(1)-0.0);
+
+%         mu = 1.5; beta = 1;
+%         yp(2) = -mu*y(2)+beta*y(1)*y(2);     %SIR
+%         yp(1) = -beta*y(1)*y(2);
+
+%      yp(1) = w*(y(2)-y(3));
+%      yp(2) = w*(y(3)-y(1));
+%      yp(3) = w*(y(1)-y(2));
+
+%     yp(1) = y(1)*(1 - 2*y(1) + y(2));
+%     yp(2) = y(2)*(1 + y(1) - y(2));
+
+        %yd = [yp(1);yp(2);yp(3)];
+        yd = [yp(1);yp(2)];
+
+    end     % end f5
+
+end
+
diff --git a/rksimpleVirus2.m b/rksimpleVirus2.m
new file mode 100644
index 0000000..ca137ca
--- /dev/null
+++ b/rksimpleVirus2.m
@@ -0,0 +1,47 @@
+
+function rksimpleVirus2
+
+r = 2.4;
+a = 2;
+b = 0.1;
+c = 1;
+
+
+x0 = r/a;
+
+
+
+xstar = r/a
+vstar = b*r/(a*c)
+
+J = [-b c; -a*vstar r-a*xstar]
+eig(J)
+Delta = det(J)
+tau = trace(J)
+
+
+
+y0 = [0.3509 1];
+tspan = [0 200];
+
+[t,y] = ode45(@f5,tspan,y0);
+
+figure(1)
+plot(t,y(:,1),t,y(:,2))
+
+figure(2)
+plot(y(:,1),y(:,2))
+
+    function yd = f5(t,y)
+
+        yp(1) = c*y(2) - b*y(1);
+        yp(2) = y(2).*(r - a*y(1));
+        
+
+        
+        yd = [yp(1);yp(2)];
+
+    end     % end f5
+
+end
+
diff --git a/rossler.m b/rossler.m
new file mode 100644
index 0000000..3ccd8b7
--- /dev/null
+++ b/rossler.m
@@ -0,0 +1,93 @@
+% function rossler
+
+function rossler
+
+a = 0.15;        % 0.15 (defined phase)     0.3 (chaotic phase)
+b = 0.4;         % 0.4
+c = 8;       % 6 8.0 
+
+y0 = [1 1 0];
+tspan = [1 30];
+[t,y] = ode45(@f5,tspan,y0);
+[sz dum] = size(y);
+
+y0 =[y(sz,1) y(sz,2) y(sz,3)];
+tspan = [1 200];
+options = odeset('OutputFcn',@odephas3);
+[t,y] = ode45(@f5,tspan,y0,options);
+
+figure(2)
+plot(t,y(:,1),t,y(:,2),t,y(:,3))
+
+figure(20)
+colormap(jet)
+plot3(y(:,1),y(:,2),y(:,3),'k')
+set(gcf,'Color','white')
+grid on
+
+
+figure(3)
+plot(y(:,1),y(:,3))
+title('x-z')
+
+figure(4)
+plot(y(:,2),y(:,3))
+title('y-z')
+
+figure(5)
+plot(y(:,1),y(:,2))
+title('x-y')
+
+PowS = 1;
+% Power Spectrum
+if PowS == 1        % Turn this off by setting PowS = 0
+    [st dum] = size(t);
+    tshort = t(400:st,1);
+    yshort = y(400:st,1);
+    %yshort = u(400:st,1);
+    [st dum] = size(tshort);
+    dt = (max(tshort)-min(tshort))/st;
+    df = 1/st/dt;
+    fmax = 1/dt;
+    freq = df:df:fmax;
+
+    yf = yshort(:,1);
+    F = fft(yf);
+    Pow = F.*conj(F);
+
+    mid = round(st/4);
+    figure(6)
+    semilogx(freq(2:mid),Pow(2:mid))
+    title('Power Spectrum')
+    xlabel('Frequency (Hz)')
+    
+
+    [Pmax I] = max(Pow);
+    f0 = freq(I)
+    T0 = 1/f0
+%     T0nat = 2*pi/sqrt(w02)
+end
+
+
+
+d = y(:,1);
+e = y(:,2);
+f = y(:,3);
+Printfile4('lout',t',d',e',f');
+
+%keyboard
+
+
+    function yd = f5(t,y)
+
+        yp(1) = -(y(2) + y(3));
+        yp(2) = y(1) + a* y(2);
+        yp(3) = b + y(3)*(y(1) - c);
+
+        yd = [yp(1);yp(2);yp(3)];
+
+    end     % end f5
+
+
+end % end ltest
+
diff --git a/sinecircle1.m b/sinecircle1.m
new file mode 100644
index 0000000..ab96eee
--- /dev/null
+++ b/sinecircle1.m
@@ -0,0 +1,77 @@
+%sinecircle1.m
+
+clear
+
+%omega = 1.618033988749895;   %exp(1)/2     pi/2     golden
+%omega = 3/2;
+omega = 0;
+%omega = 1;
+%omega = sqrt(2);
+%omega = 21/13;
+%omega = 13/8;
+%omega = 8/5;
+%omega = 5/3;
+%omega = 2/3;
+
+testcase = 2;           % 1 = iterate at a single g     2= scan g
+
+if testcase == 1
+    gval = 80.;         % 43 76
+    glooplo = gval;
+    gloophi = gval;
+elseif testcase == 2   % scan g
+    figure(3)
+    close
+    glooplo = 0;
+    gloophi = 2;
+end
+
+
+delg = 0.005;
+for gloop = glooplo:delg:gloophi
+
+    g =gloop;
+
+    if testcase == 1
+        disp(strcat('g = ',num2str(g)))
+    end
+
+    thet0 = rand;
+    thet1last = thet0;
+    N = 100;
+    for loop = 1:N
+
+        thet1emp = mod(thet1last + omega + g*sin(2*pi*thet1last),1);
+        thet1last = mod(thet1emp,1);
+
+    end
+
+    for loop = 1:N
+
+        thet1(loop) = mod(thet1last + omega + g*sin(2*pi*thet1last),1);
+        thet1last = mod(thet1(loop),1);
+
+    end
+
+    if testcase == 1
+        figure(1)
+        plot(thet1,thet2,'o')
+        axis([0 1 0 1])
+
+        figure(2)
+        plot(thet1,(thet2-thet1),'o')
+        axis([0 1 0 1])
+    elseif testcase == 2
+        x = g*ones(1,N);
+        y = thet1;
+        hold on
+        figure(3)
+        plot(x,y,'.k')
+        axis([0 gloophi 0 1])
+        hold off
+        set(gcf,'Color','White')
+    end
+
+end
+
+
diff --git a/zerodiag.m b/zerodiag.m
new file mode 100644
index 0000000..765ac5b
--- /dev/null
+++ b/zerodiag.m
@@ -0,0 +1,24 @@
+% [B, D] = zerodiag(A,altmin,varargin)
+% Sets diagonal terms to zero
+% D is the diagonal values
+
+function [B, D] = zerodiag(A,altmin,varargin)
+
+[sy,sx] = size(A);
+B = A;
+if sy~=sx
+    disp('in zerodiag not square');
+    B = A;
+else
+    if nargin == 1
+        for loop = 1:sx
+            B(loop,loop) = 0;
+            D(loop) = A(loop,loop);
+        end
+    else
+        for loop = 1:sx
+            B(loop,loop) = altmin;
+            D(loop) = A(loop,loop);
+        end
+    end
+end
\ No newline at end of file