function ode2ndOrder
clear all
tr=[0 15]; %seconds
initv=[600 0]; %start 600 m high
[t,y]=ode45(@rkfalling, tr, initv);
plot(t,y(:,1))
ylabel('x (m)')
xlabel('time(s)')
figure
plot(t,y(:,2))
ylabel('velocity (m/s)')
xlabel('time(s)')

function r=rkfalling(t,w)
mass=80;
k=4/15/mass;
g=9.81;
r=zeros(2,1);
r(1)=w(2);
r(2)= k*w(2)^2-g;