clear
global Omegar Omegam OmegaL Omega0 alpha

irn=4;
alpha=1;
tmax=100;

switch irn
    case 1   % flat matter only
        Omega0=1;
        Omegar=0;
        OmegaL=0;
        Omegam=Omega0-(Omegar+OmegaL);
    case 2   % neg. curved matter only
        Omega0=0.8;
        Omegar=0;
        OmegaL=0;
        Omegam=Omega0-(Omegar+OmegaL);
    case 3   % pos. curved matter only
        Omega0=1.1;
        Omegar=0;
        OmegaL=0;
        Omegam=Omega0-(Omegar+OmegaL);
    case 4  % Benchmark model
        Omega0=1;
        Omegar=8.4e-5;
        Omegam=0.30;
        OmegaL=Omega0-(Omegar+Omegam);
        tmax=10;
    case 5   % neg. curved radiation only
        Omega0=0.8;
        Omegam=0;
        OmegaL=0;
        Omegar=Omega0-(Omegar+OmegaL);
    case 6   % pos. curved radiation only
        Omega0=0;                     % 1.1;
        Omegam=0;
        OmegaL=0;
        Omegar=Omega0-(Omegar+OmegaL);
        %text(60,400,['\Omega_0 =  ',num2str(0.8)])
        %text(70,200,['\Omega_0 =  ',num2str(0)])
        %text(80,50,['\Omega_0 =  ',num2str(1.1)])
    otherwise
end    

if irn<5
    
    f=@(x)  Omegar./x.^2+Omegam./x+OmegaL*x.^2+1-Omega0;
    amax=fzero(f,10)

    figure(1)
    tspan=0:0.1:tmax;
    a0=.1;
    [t,a]=ode45(@univf,tspan,a0);

    if (amax>1 & isnan(amax)==0) 
        ind=(find(a>amax));
        tmax=t( ind(1)-1);
        tspan=0:4:tmax;
        a0=.1;
        [t,a]=ode45(@univf,tspan,a0);
    end

    plot(t,real(a),'ok','MarkerSize',2)
    xlabel('t')
    ylabel('a')

    switch irn
        case 1
            hold on
            y=(3*t/2).^(2/3);
            plot(t,y,'b')
            hold off
        case 2
            hold on
            eta=0:.1:4;
            a3 = Omega0/(1-Omega0)/2*(cosh(eta)-1);
            t3 = Omega0/(1-Omega0)^1.5/2*(sinh(eta)-eta);
            plot(t3,a3,'k')
            axis([0,max(t3),0,max(a3)])
            hold off
        case 3
            if amax>1
                a0=amax-.001;
                tspan=tmax:4:1.8*tmax;
                alpha=-1;
                [t,a]=ode45(@univf,tspan,a0);
                hold on
                plot(t,real(a),'ok','MarkerSize',2)
                hold off
            end
            hold on
            phi=0:.1:2*pi;
            a3 = Omega0/(Omega0-1)/2*(1-cos(phi));
            t3 = Omega0/(Omega0-1)^1.5/2*(phi-sin(phi));
            plot(t3,a3,'k')
            hold off
        otherwise
    end
    
elseif irn<7
    f=@(x)  Omegar./x.^2+Omegam./x+OmegaL*x.^2+1-Omega0;
    figure(1)
    tspan=0:1:tmax;
    a0=.1;
    [t,a]=ode45(@univf,tspan,a0);
    plot(t,real(a),'ok','MarkerSize',2)
    xlabel('t')
    ylabel('a')
else
    Omegar=1e-05;
    Omegam=1-Omegar;
    OmegaL=0;
    f=@(x)  Omegar./x.^2+Omegam./x;
    
    arm=Omegar/Omegam;
    amax=fzero(f,10);
    tmax=1;

    figure(1)
    tspan=0:.01:tmax;
    a0=arm/100;
    [t,a]=ode45(@univf,tspan,a0);
    
    plot(t,real(a),'ok','MarkerSize',2)
    xlabel('t')
    ylabel('a')
    hold on
    a4=0:.1:1.3;
    t4 = 4*arm^2*(1-(1-a4/(2*arm)).*(1+a4/arm).^(1/2))/(3*sqrt(Omegar));
    plot(t4,a4,'b')
    hold off
end
    
