global c q m;

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


dt = 0.004;

t(1)=0;
y(1,1)=-60;
y(1,2)=0;
z=forma(y(1,1));
eu(1)=9.8*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);
av(1)=0.;

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

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

z=forma(y(i+1,1));
eu(i+1)=9.8*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);
av(i+1)=(zz(i+1)-zz(i))/dt;

endfor