Matlab-Friedmann Equations, help?

[birdhen]

I have been asked to create a programme on Matlab to solve the Friedmann equations. However being quite computer illiterate I am not really sure where to begin. I have been reading some books on Matlab, but still not sure how to go about this. Could anyone offer any tips?

 

[fab13]

here is a numerical solving program with MATLAB "ode" solvers,
most useful for ΛCDM models where there is not analytic solution :

function friedmann

H0=71000/(3*10^(22))*(3600*24*365*10^(9));

t_begin(1)=13.7;
t_begin(2)=13.7;
t_final(1)=0;
t_final(2)=40;
for j=1:2
    
Omega1_m=0.3;
Omega1_red=0;
Omega1_k=1-Omega1_m-Omega1_red;
Omega2_m=0.3;
Omega2_red=0.7;
Omega2_k=1-Omega2_m-Omega2_red;
Omega3_m=5;
Omega3_red=0;
Omega3_k=1-Omega3_m-Omega3_red;
Omega4_m=1;
Omega4_red=0;
Omega4_k=1-Omega4_m-Omega4_red;

a1_0=1;
a1_prim_0=H0*(Omega1_m/a1_0+Omega1_k)^(1/2);
init1=[ a1_0 ; a1_prim_0 ];
a2_0=1;
a2_prim_0=H0*(Omega2_m/a2_0+Omega2_k+Omega2_red*a2_0)^(1/2);
init2=[ a2_0 ; a2_prim_0 ];
a3_0=1;
a3_prim_0=H0*(Omega3_m/a3_0+Omega3_k)^(1/2);
init3=[ a3_0 ; a3_prim_0 ];
a4_0=1;
a4_prim_0=H0*(Omega4_m/a4_0+Omega4_k)^(1/2);
init4=[ a4_0 ; a4_prim_0 ];
options = odeset('Events',@events,'RelTol',10^(-10),'AbsTol',10^(-10));


[t1,y1]=ode45(@(t1,y1) syseq(t1,y1,Omega1_m,Omega1_red,Omega1_k),[t_begin(j) t_final(j)],init1,options);
[t2,y2]=ode45(@(t2,y2) syseq(t2,y2,Omega2_m,Omega2_red,Omega2_k),[t_begin(j) t_final(j)],init2,options);
[t3,y3]=ode45(@(t3,y3) syseq(t3,y3,Omega3_m,Omega3_red,Omega3_k),[t_begin(j) t_final(j)],init3,options);
[t4,y4]=ode45(@(t4,y4) syseq(t4,y4,Omega4_m,Omega4_red,Omega4_k),[t_begin(j) t_final(j)],init4,options);

figure(1);
hold on;
plot(t1,y1(:,1),'b',t2,y2(:,1),'r',t3,y3(:,1),'k',t4,y4(:,1),'g'); 

xlabel('Time in Gyr');
ylabel('Relative size of the universe');
ylim([ 0 4]);

end
end

function dydt = syseq(t,y,Omega_m,Omega_red,Omega_k) 

H0=71000/(3*10^(22))*(3600*24*365*10^(9));

dydt=[ y(2) ;                        
 (-(H0^(2)/2)*(Omega_m/y(1)^(3)-2*Omega_red))*Omega_m*(1/(H0^(2))*y(2)^2-Omega_k-Omega_red*y(1)^(2))^(-1);                
       ];
 
   
end


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);     % Detect height = 0
isterminal = 1;   % Stop the integration
direction = -1;   % Negative direction only
end

Regards