diff --git a/SIRSMC.m b/SIRSMC.m
new file mode 100644
index 0000000..9b052fa
--- /dev/null
+++ b/SIRSMC.m
@@ -0,0 +1,257 @@
+% SIRSMC.m
+% converted from diffusionmat on 1/4/16
+% Susecpt-Infect-Suscept  plus an innoculation of highest-degree node
+% Monte Carlo to get good statistics
+
+clear
+close all
+format compact
+
+T1 = 100;
+T2 = 300;
+
+showfig = 0;
+
+Innoc = 1;  % Turn innoculation on
+
+Nensemble = 100;
+for eloop = 1:Nensemble
+    eloop
+    
+    % beta is infection rate
+    % mu is recovery rate
+    
+    %N = 50;  m = 2;    beta = 0.2;  mu = 0.6;               node = makeSF(N,m);
+    %N = 50;  m = 2;    beta = 0.2;  mu = 0.4;   p = 0.1;   node = makeSW(N,m,p);
+    N = 50;  p = 0.06;  beta = 0.2;  mu = 0.4;               node = makeER(N,p);
+    
+    [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(' ')
+    disp(strcat('fac =',num2str(mu/avgdegree/beta)))
+    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
+    a = zeros(N,1);
+    [Y,I] = max(degree);
+    a(I) = 1;   % Infect the high-degree node
+    
+    if showfig == 1
+        fh2 = figure(2);
+        drawnet(node)
+    end
+    
+    % The discrete-time approach
+    
+    c0 = a;
+    dt = 1;
+    
+    R = eye(N,N);
+    Ct = zeros(1,N);
+    for nloop = 1:N
+        deg(nloop) = node(nloop).numlink;
+    end
+    
+    tic
+    c = c0;
+    fh3 = figure(3);
+    eps = 1e-6;
+    Con = zeros(T2,N);
+    in = zeros(1,T2);
+    flag = 0; timeloop = 0;
+    while (flag == 0)&&(timeloop <= T1)
+        timeloop = timeloop + 1;
+        
+        %M = eye(N,N) - beta*Lap*dt.*randbin2(N,N,1-beta).*(ones(N,N)-eye(N,N));
+        M = eye(N,N) + beta*A*dt.*randbin2(N,N,1-beta).*(ones(N,N)-eye(N,N));
+        
+        ctmp = M*c;
+        
+        ctmp2 = ceil(ctmp);
+        
+        
+        c = maskbilevel(ctmp2,0,1,0,1);
+        
+        
+        ctmp3 = R*c;
+        ctmp4 = floor(ctmp3);
+        c = maskbilevel(ctmp4,0,1,0,1);
+        
+        
+        Con(timeloop,:) = c';
+        
+        Pop = sum(Con(timeloop,:));
+        if Pop == 0
+            flag = 1;
+        end
+        
+        
+        for nodeloop = 1:N
+            node(nodeloop).value = c(nodeloop);
+            Rtmp(nodeloop,nodeloop) = c(nodeloop);
+        end
+        
+        if showfig == 1
+            drawnet(node,2)
+            pause(0.01)
+        end
+        
+        for nloop = 1:N
+            if node(nloop).value == 1
+                Ct(nloop) = Ct(nloop) + 1/T1;
+            end
+        end
+        
+        R = eye(N,N) - (Rtmp.*mu*dt.*eye(N,N).*randbin2(N,N,1-mu));
+        
+        %     if timeloop > T1-2
+        %         keyboard
+        %     end
+        
+    end
+    
+    if Innoc == 1
+        % Innocluate
+        % Remove the highest-degree node
+        displine('avgdegree = ',avgdegree)
+        disp(strcat('fac =',num2str(mu/avgdegree/beta)))
+        node = subnode(I,node);
+        %snode = removenode(I,node);
+        [A,degree,Lap] = adjacency(node);
+        
+        %keyboard
+        
+        [N,e,avgdegree,maxdegree,mindegree,numclus,meanclus,Lmax,L2,LmaxL2] = clusterstats(node);
+        displine('avgdegree = ',avgdegree)
+        disp(strcat('fac =',num2str(mu/avgdegree/beta)))
+    end
+        
+        
+        while (flag == 0)&&(timeloop <= T2)
+            timeloop = timeloop + 1;
+            
+            %M = eye(N,N) - beta*Lap*dt.*randbin2(N,N,1-beta).*(ones(N,N)-eye(N,N));
+            M = eye(N,N) + beta*A*dt.*randbin2(N,N,1-beta).*(ones(N,N)-eye(N,N));
+            
+            ctmp = M*c;
+            
+            ctmp2 = ceil(ctmp);
+            c = maskbilevel(ctmp2,0,1.01,0,1);
+            
+            ctmp3 = R*c;
+            ctmp4 = floor(ctmp3);
+            c = maskbilevel(ctmp4,-0.01,1.01,0,1);
+            
+            
+            Con(timeloop,:) = c';
+            Pop = sum(Con(timeloop,:));
+            if Pop == 0
+                flag = 1;
+            end
+            
+            
+            for nodeloop = 1:N
+                node(nodeloop).value = c(nodeloop);
+                Rtmp(nodeloop,nodeloop) = c(nodeloop);
+            end
+            
+            if showfig == 1
+                drawnet(node,2)
+                pause(0.01)
+            end
+            
+            R = eye(N,N) - (Rtmp.*mu*dt.*eye(N,N).*randbin2(N,N,1-mu));
+            
+            %     if timeloop > T1-2
+            %         keyboard
+            %     end
+            
+        end
+        toc
+        
+        
+        
+        x = 0:T2-1;
+        h = colormap(jet);
+        figure(4)
+        for tloop = 1:T2
+            Ssum = 0;
+            for nodeloop = 1:N
+                Ssum =  Ssum + Con(tloop,nodeloop);
+            end
+            In(tloop) = Ssum;
+        end
+        
+        figure(4)
+        plot(In)
+        title('Infected')
+        
+        
+        mn1 = mean(In(10:T1));
+        mn2 = mean(In(T1+10:T2));
+        
+        displine('del pop = ',mn2-mn1)
+        displine('rel del pop = ',(mn2-mn1)/mn1)
+        
+        figure(5)
+        plot(deg,Ct,'o')
+        xlabel('degree')
+        ylabel('avg infection')
+        
+        
+%     else  % non-Innoc case
+%         
+%         
+%         x = 0:T1-1;
+%         h = colormap(jet);
+%         figure(4)
+%         for tloop = 1:T1
+%             Ssum = 0;
+%             for nodeloop = 1:N
+%                 Ssum =  Ssum + Con(tloop,nodeloop);
+%             end
+%             In(tloop) = Ssum;
+%         end
+%         
+%         figure(4)
+%         plot(In)
+%         title('Infected')
+%         
+%     end
+    
+    Infection(eloop,:) = In;
+    
+end     % end eloop
+
+figure(6)
+imagesc(Infection)
+
+Y = mean(Infection);
+
+figure(7)
+plot(Y)
+
+