diff --git a/Heiles2.m b/Heiles2.m
new file mode 100644
index 0000000..5a72822
--- /dev/null
+++ b/Heiles2.m
@@ -0,0 +1,150 @@
+function Heiles2
+
+c(1) = 'r';
+c(2) = 'b';
+c(3) = 'g';
+c(4) = 'k';
+c(5) = 'c';
+c(6) = 'm';
+c(7) = 'y';
+
+h = colormap(lines);
+
+% Either E = 1 and epsE = 0.36
+% or E = 0.13 and Epse = 1
+
+
+E = .1;        % 1/12     .19
+
+epsE = 1;         % 1.0
+
+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(jet)
+imagesc(V)
+hold on
+contour(V,30,'LineColor','k')
+hold off
+caxis([0 1])
+colorbar
+
+
+figure(5)
+clf
+hold on
+repnum = 64;                % 32
+mulnum = 64/repnum;
+for reploop = 1:repnum
+    
+    clear yvar py
+    
+%     py = 0;
+%     yp1 = -0.0945;
+%     px = 2*(rand-0.5)*sqrt(2*(E-py^2/2-yp1^2/2));
+%     xp1 = sign(rand-0.5)*real(sqrt(2*(E-px^2/2-yp1^2/2-py^2/2)));
+
+%Random split E between px and py
+    px = 2*(rand-0.499)*pmax;
+    py = sign(rand-0.5)*real(sqrt(2*(E-px^2/2)));
+    xp1 = 0;
+    yp1 = 0;
+
+    Etest = px^2/2 + py^2/2 + xp1^2/2 + yp1^2/2
+        
+    y0 = [xp1 yp1 px py];
+    
+    tspan = [1 2000];            % 500
+
+    figure(1)
+    %options = odeset('OutputFcn',@odeplot,'RelTol',1e-7);
+    options = odeset('OutputFcn',@odephas3,'RelTol',1e-7);
+    % options = odeset('OutputFcn',@odeplot);
+
+    [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;
+                py(cnt) = y4(loop);
+                yvar(cnt) = y2(loop);
+                last = y1(loop);
+            else
+                last = y1(loop);
+            end
+        end
+
+        figure(5)
+        plot(yvar,py,'o','MarkerSize',4,'Color',h(floor(0.9*mulnum*reploop)+1,:),'MarkerFaceColor',h(floor(0.9*mulnum*reploop)+1,:))
+        xlabel('y')
+        ylabel('py')
+
+    end
+
+end     % end reploop
+
+figure(5)
+hold off
+
+
+% Model function
+    function dy = f1(t,y)
+
+        B = 0.000;          %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/Lozi2.m b/Lozi2.m
new file mode 100644
index 0000000..114d9d6
--- /dev/null
+++ b/Lozi2.m
@@ -0,0 +1,46 @@
+% Lozi.m
+
+clear
+
+B = -1.00;
+C = 0.5;
+
+%h = colormap(hsv);
+h = colormap(jet);
+
+f1 = figure(1);
+f1.Position = [382 147 976 808];
+dum = set(f1);
+
+for itloop = 1:50    %200
+    itloop
+    rn = ceil(255*rand);
+    
+    xlast = randn;
+    ylast = randn;
+    
+%     for transloop = 1:100   % Transient loop not needed for Hamiltonian case
+%         xnew = 1 + ylast - C*abs(xlast);
+%         ynew = B*xlast;
+%         xlast = xnew;
+%         ylast = ynew;
+%     end
+    
+    figure(1)
+    %axis([-160.5 160.5 -160.5 160.5])
+    axis([-2 2 -2 2])
+    %clf
+    hold on
+    for loop = 1:500
+        xnew = 1 + ylast - C*abs(xlast);
+        ynew = B*xlast;
+        xlast = xnew;
+        ylast = ynew;
+        plot(real(xnew),real(ynew),'o','MarkerSize',3,'LineWidth',0.5,'Color',h(rn,:))
+    end
+    %plot(real(xnew),real(ynew),'o','MarkerSize',2,'Color','r')
+    pause(0.001)
+    
+    %hold off
+    
+end
\ No newline at end of file