function friedmann

clear;
clf;
%%%% Actual Hubble value (m.Gyr^-1)
H0=67400/(3.09*10^(22))*(3600*24*365*10^(9));
%%%% Starting and ending time
t_begin(1) = 0;
t_final(1) = -15;
t_begin(2) = 0;
t_final(2) = 30;
%%%% Initial Omega parameters    
Omega1_m=0.3;
Omega1_r=10^-4;
Omega1_l=0;
Omega1_k=1-Omega1_m-Omega1_r-Omega1_l;
Omega2_m=0.3;
Omega2_r=10^-4;
Omega2_l=0.7;
Omega2_k=1-Omega2_m-Omega2_r-Omega2_l;
Omega3_m=5;
Omega3_r=10^-4;
Omega3_l=0;
Omega3_k=1-Omega3_m-Omega3_r-Omega3_l;
Omega4_m=1;
Omega4_r=10^-4;
Omega4_l=0;
Omega4_k=1-Omega4_m-Omega4_r-Omega4_l;
%%%% Initial scale factor and derivated scale factor
a1_0=1;
a1_prim_0=H0*(Omega1_m/a1_0+Omega1_r/(a1_0^2)+Omega1_k+Omega1_l*(a1_0^2))^(1/2);
init1=[ a1_0 ; a1_prim_0 ];
a2_0=1;
a2_prim_0=H0*(Omega2_m/a2_0+Omega2_r/(a2_0^2)+Omega2_k+Omega2_l*(a2_0^2))^(1/2);
init2=[ a2_0 ; a2_prim_0 ];
a3_0=1;
a3_prim_0=H0*(Omega3_m/a3_0+Omega3_r/(a3_0^2)+Omega3_k+Omega3_l*(a3_0^2))^(1/2);
init3=[ a3_0 ; a3_prim_0 ];
a4_0=1;
a4_prim_0=H0*(Omega4_m/a4_0+Omega4_r/(a4_0^2)+Omega4_k+Omega4_l*(a4_0^2))^(1/2);
init4=[ a4_0 ; a4_prim_0 ];
set(groot,'defaultLineLineWidth',2.0)
%%%% Solving in two parts : backwards and forwards
for i=1:2
 %%%% options for solver
 options = odeset('Events',@events,'RelTol',10^(-10),'AbsTol',10^(-10));
 %%%% Differential systems solving
 [t1,y1]=ode45(@(t1,y1) syseq(t1,y1,Omega1_m,Omega1_r,Omega1_l,Omega1_k),[t_begin(i) t_final(i)],init1,options);
 [t2,y2]=ode45(@(t2,y2) syseq(t2,y2,Omega2_m,Omega2_r,Omega2_l,Omega2_k),[t_begin(i) t_final(i)],init2,options);
 [t3,y3]=ode45(@(t3,y3) syseq(t3,y3,Omega3_m,Omega3_r,Omega3_l,Omega3_k),[t_begin(i) t_final(i)],init3,options);
 [t4,y4]=ode45(@(t4,y4) syseq(t4,y4,Omega4_m,Omega4_r,Omega4_l,Omega4_k),[t_begin(i) t_final(i)],init4,options);
 figure(1);
 hold on;
 %%%% Plot solutions
 plot(t1,y1(:,1),'b',t2,y2(:,1),'r',t3,y3(:,1),'k',t4,y4(:,1),'m');
 box on;
 grid on;
 xlabel('Cosmic Time in Gyr');
 ylabel('Relative size of the universe');
 ylim([ 0 4]);
 %%%% Set label for each curve 
 posX = -13;
 posX1 = posX;
 posX2 = posX;
 posX3 = posX;
 posX4 = posX;
 posY1 = 3.6;
 posY2 = 3.2;
 posY3 = 2.8;
 posY4 = 2.4;
 strMax = ['$$\Omega^{0}_{m} = 0.3 , \Omega^{0}_{\Lambda} = 0 , \Omega^{0}_{r} = 10^{-4} , k = -1$$']; 
 text(posX1, posY1, strMax, 'FontSize', 20, 'Color', 'blue', 'VerticalAlignment', 'top', 'Interpreter', 'latex');
 strMax = ['$$\Omega^{0}_{m} = 0.3 , \Omega^{0}_{\Lambda} = 0.7 , \Omega^{0}_{r} = 10^{-4} , \Omega^{0}_{k} \simeq 0$$'];
 text(posX2, posY2, strMax, 'FontSize', 20, 'Color', 'red', 'VerticalAlignment', 'top', 'Interpreter', 'latex');
 strMax = ['$$\Omega^{0}_{m} = 5 , \Omega^{0}_{\Lambda} = 0 , \Omega^{0}_{r} = 10^{-4} , k = 1$$'];
 text(posX3, posY3, strMax, 'FontSize', 20, 'Color', 'black', 'VerticalAlignment', 'top', 'Interpreter', 'latex');
 strMax = ['$$\Omega^{0}_{m} = 1 , \Omega^{0}_{\Lambda} = 0 , \Omega^{0}_{r} = 10^{-4} , k = 1$$'];
 text(posX4, posY4, strMax, 'FontSize', 20, 'Color', 'magenta', 'VerticalAlignment', 'top', 'Interpreter', 'latex');
 title('Relative size of the universe vs Cosmic Time');
end
end
%%%% Differential system function
function dydt = syseq(t,y,Omega_m,Omega_r,Omega_l,Omega_k) 
%%%% Actual Hubble value
H0=67400/(3.09*10^(22))*(3600*24*365*10^(9));
dydt=[ y(2) ;
      (-(H0^(2)/2)*(Omega_m/y(1)^(3)+2*Omega_r/y(1)^(4)-2*Omega_l))*Omega_m*(1/(H0^(2))*y(2)^(2)-...
      Omega_r/y(1)^(2)-Omega_l*y(1)^(2)-Omega_k)^(-1);];
end
%%%% Integration stopping function
function [value,isterminal,direction] = events(t,y)
% Locate the time when height passes through zero in a
% decreasing direction and stop integration.
value = [y(1),4-y(1)]; % Detect height = 0
isterminal = [1,1];    % Stop the integration
direction = [-1,1];    % Negative and positive directions
end