clc, clear all, close all
	y0=[1;-1.57];xini=1e-15;xfin=60;
	[x,y] = ode113(@vdp1,[xini xfin],y0);
	de=1e-4;while y(length(y),2) > 0,y0(2,1)=y0(2,1)-de;
	[x,y] = ode113(@vdp1,[xini xfin],y0);
	if y(length(y),2) < 0,y0(2,1)=y0(2,1)+de;
	[x,y] = ode113(@vdp1,[xini xfin],y0);
	de=de/2;elseif y(length(y),2) < 0.0001,
	y(length(y),2),plot(x,y(:,1),'.k','MarkerSize',15)
	xlabel('x','FontSize',18),ylabel('\chi(x)','FontSize',18)
	set(gca,'FontSize',18),axis( [ 0 25 0 1 ] );
	save tempX x,yy=y(:,1);save tempYeYl yy,return,end,end

	function dydt = vdp1(x,y)
	dydt = [y(2); y(1).^(3/2)./(x.^(0.5))];


	clc,clear all,close all
	load tempX,load tempYeYl  
	Z=2; %% Definindo número atômico
	cf=(3/10)*(3*pi^2)^(2/3);
	alpha=((4*pi)^(2/3))*(3/(5*cf))*Z^(1/3);
	ChiX=yy./x;Rho=(32/(9*pi^3))*(ChiX).^(3/2)*Z^2;
	norm=trapz(x,Rho.*x.^2)*4*pi/alpha^3 %norm deve ser igual a Z
	kernelcin=4*pi*Rho.^(5/3).*x.^2/alpha^3;
	Ecin=cf*trapz(x,kernelcin) %%  Energia cinética
	kernelVen = -Z*4*pi*Rho.*x/alpha.^2;%% E elétron-núcleo
	Ven=trapz(x,kernelVen),for i=1:size(Rho,1);for j=1:size(Rho,1);
	Rho12(i,j)=Rho(i)*Rho(j);end,end,x1=x;x2=x';
	for i=1:size(x1,1);for j=1:size(x2,2);
	X12(i,j)=(x1(i)+x2(j)-abs(x1(i)-x2(j)))*x1(i)*x2(j);
	end,end,kernelVee=(4*(pi)^2/(alpha^5))*Rho12.*X12;
	Vee=trapz(x2,trapz(x1,kernelVee,1),2) %% E elétron-elétron
	Etot=Ecin+Ven+Vee 	%% Energia total