global c q m;

c=4;
q=-1;
m=1;


dt = 0.001;

t(1)=0;
y(1,1)=-36;
y(1,2)=0;
z=forma(y(1,1));
eu(1)=z(1);
ek(1)=(1+z(2)^2)*y(1,2)^2/2;
e(1)=eu(1)+ek(1);
zz(1)=z(1);
dh(1)=z(2);

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

for j=1:5
	yy=yy+dt*fmrusa(tt,yy);
	tt=tt+dt;
endfor
t(i+1)=tt;
y(i+1,:)=yy;

z=forma(y(i+1,1));
eu(i+1)=z(1);
ek(i+1)=(1+z(2)^2)*y(i+1,2)^2/2;
e(i+1)=eu(i+1)+ek(i+1);
zz(i+1)=z(1);
dh(i+1)=z(2);

endfor