function [t,y] = gravi_edo(peri,afel)
	GM = 1.327124e20;
	h = sqrt(2*GM/(peri^-1+afel^-1));
	% periodo con la tercera ley de Kepler
	period = sqrt((peri+afel)^3*pi^2/2/GM);
	[t,y] = ode45( @eqn, [0,period], [peri,0,0]);
	% EDO de la p.2 del resumen
	function yp = eqn(t,y)
	%	yp = [-(y(2)); y(1)];
		rp = y(2);
		sp =(h^2/y(1)-GM)/y(1)^2;
		thp = h/y(1)^2;
		yp = [rp; sp; thp];
	end
end
