global c q m;

c=1;
q=-1;
m=1;
centro=[0;0];

dt = 0.005;

t(1)=0;
y(1,1)=-4;
y(1,2)=-5;
y(1,3)=.4;
y(1,4)=0;
r=sqrt(y(1,1)^2+y(1,2)^2);
v=sqrt(y(1,3)^2+y(1,4)^2);
e(1)=c*q/r+m*v*v/2;
ek(1)=m*v*v/2;
eu(1)=c*q/r;

plot(y(1,1),y(1,2),"box","markersize",15,centro(1),centro(2),"go","markersize",20);
axis ([-20 10 -15 15]);
hold on;

for i=1:5000
yy=y(i,:);
tt=t(i);

for j=1:40
	yy=yy+dt*felet(tt,yy);
	tt=tt+dt;
#plot(yy(1),yy(2),"ro","markersize",10);
endfor
t(i+1)=tt;
y(i+1,:)=yy;
r=sqrt(y(i+1,1)^2+y(i+1,2)^2);
v=sqrt(y(i+1,3)^2+y(i+1,4)^2);
e(i+1)=c*q/r+m*v*v/2;
ek(i+1)=m*v*v/2;
eu(i+1)=c*q/r;
#plot(y(i+1,1),y(i+1,2),"ro","markersize",8);
#pause;
endfor