clear
global Omega
figure(1)
hold on

for Omega=0.8:.1:1.2;
    if Omega<1
        amax=50;
        tmax=100;
    elseif Omega==1
        amax=50;
        tmax=100;
    else
        amax=Omega/(Omega-1);
        tmax=Omega/(Omega-1)^1.5/2*pi;
    end


    tspan=0:4:tmax;
    a0=.1;
    [t,a]=ode45(@cosmosf,tspan,a0);
    plot(t,a,'ok','MarkerSize',2)
    hold on
    if Omega>1
         tspan=tmax+.0001:4:2*tmax;
         a0=amax-.0001;
         [t2,a2]=ode45(@cosmosf2,tspan,a0);
         plot(t2,a2,'ok','MarkerSize',2)
    end

   % N=10;
   % da=amax/N;
   % F=@(a) 1./sqrt(Omega./a+(1-Omega));
   % for i=1:N
   %     a2(i)=i*da;
   %     Ht(i)=quadgk(F,0,a2(i));
   % end
   % plot(Ht,a2,'ok')

    if Omega<1
        eta=0:.1:4;
        a3 = Omega/(1-Omega)/2*(cosh(eta)-1);
        t3 = Omega/(1-Omega)^1.5/2*(sinh(eta)-eta);
        plot(t3,a3,'k')
        axis([0,max(t3),0,max(a3)])
        xlabel('t')
        ylabel('a')
    elseif Omega==1
        t3=0:.1:1.5*tmax;
        a3=(3*t3/2).^(2/3);
        plot(t3,a3,'k')
    else    
        phi=0:.1:2*pi;
        a3 = Omega/(Omega-1)/2*(1-cos(phi));
        t3 = Omega/(Omega-1)^1.5/2*(phi-sin(phi));
        plot(t3,a3,'k')
    end
end
hold off 
axis([0,150,0,50])
xlabel('t')
ylabel('a')
text(100,7,['\Omega_0 =  ',num2str(1.1)])
text(45,2,['\Omega_0 =  ',num2str(1.2)])
text(125,30,['\Omega_0 =  ',num2str(1)])
text(125,45,['\Omega_0 =  ',num2str(0.9)])
text(60,45,['\Omega_0 =  ',num2str(0.8)])