function odescripts
clear all
timerange=[0 15]; %seconds
initialvelocity=0; %meters/second
[t,y]=ode45(@f,timerange, initialvelocity);
plot(t,y)
ylabel('velocity (m/s)')
xlabel('time(s)')
disp('Press any key to continue')
pause

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)')
disp('Press any key to continue')
pause

tspan=[0 1.1];
wnot(1)=0; wnot(2)=10;
wnot(3)=0; wnot(4)=10;
[t,y]=ode45(@rkprojectile,tspan,wnot);
plot(t,y(:,1),t,y(:,3))
figure
plot(y(:,1),y(:,3))

function rk=f(t,y)
mass=80;
g=9.81;
rk=-g+4/15*y^2/mass;

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;

function r=rkprojectile(t,w)
g=9.81;
x=w(1); s=w(2); y=w(3); z=w(4);
vel=sqrt(s.^2+z.^2);
r=zeros(4,1);
r(1)=s;
r(2)=-s*vel;
r(3)=z;
r(4)=-z*vel-g;