function light_geodesics

%% Input parameters
dist=input('Distance of the galaxy (Mpc): ');
%% Hubble constant
H0=71000;
%% Light velocity
cv=3*10^(8);
%% Time begin
ctk=(2*cv)/(3*H0);
w1=ctk;
%% Initial conditions
z1=pi/2;
y1=pi/2;
x1=dist;
xp=-1;
yp=0;
zp=0.;
%% Compute initial d(ct)/ds
wp=sqrt((3*H0/(2*cv)*w1)^(4/3)*(xp^(2)+x1^(2)*yp^(2)+x1^(2)*(sin(y1))^(2)*zp^(2)));
%% Options for ode solver
options=odeset('NonNegative',3);
%% Initial conditions array
init=[ w1;wp;x1;xp;y1;yp;z1;zp];
%% Solving differential system
[s,y]=ode45(@syseq,0:10000,init,options);
%% Apply scale factor for physical distance
y(:,3)=y(:,3).*(3*H0/(2*cv)*y(:,1)).^(2/3);
%% Plotting solution
figure(1);
plot(y(:,1),y(:,3));
box on;
grid on;
xlabel('Local Time (Mpc)');
ylabel('Distance (Mpc)');
ylim([0 dist]);
title('Light ray trajectory with Matlab ODE');
%% Plotting analytical solution for radial trajectory
figure(2);
Dphi=(3*H0/(2*cv)*y(:,1)).^(2/3)*dist-3*y(:,1).^(2/3).*(y(:,1).^(1/3)-ctk^(1/3));
%% Plotting solution
plot(y(:,1),Dphi);
box on;
grid on;
xlabel('Local Time (Mpc)');
ylabel('Distance (Mpc)');
ylim([0 dist]);
title('Light ray trajectory with analytical solution');
end        
%% Definition of differential system
function dydt = syseq(s,y) 
H0=71000;
cv=3*10^(8);
 dydt=zeros(8,1);
 dydt=[ y(2) ;                        
       -((2/3)*((3*H0/(2*cv))^(4/3))*y(1)^(1/3))*(y(4)^(2)+abs(y(3))^(2)*y(6)^(2)+abs(y(3))^(2)*(sin(y(5)))^(2)*y(8)^(2)) ;
       y(4) ;
       -(4/(3*y(1)))*y(2)*y(4)+abs(y(3))*(y(6)^(2)+(sin(y(5)))^(2)*y(8)^(2));                            
       y(6) ;
       -(2/abs(y(3)))*y(4)*y(6)-(4/(3*y(1)))*y(2)*y(6)+((sin(2*y(5)))*(y(8)^(2))/2);                           
       y(8) ;
       -2*y(8)*((2/(3*y(1)))*y(2)+(cot(y(5)))*y(6)+(1/abs(y(3)))*y(4)) ];                   
end