T = linspace(0,100,1000);
x0=0;
y0=0;
vx0= 20 + 79;
vy0= 20 + 79;
t =[x0,y0,vx0,vy0];
function dt=f(t,T)
b = 0.02;
g = 9.81;
dt = zeros(4,1);
dt(1) = t(3);
dt(2) = t(4);
dt(3) = -b*t(3)*sqrt((t(3))^2+(t(4))^2);
dt(4) = -b*t(4)*sqrt((t(3))^2 +(t(4))^2) - g;
endfunction
sol = lsode("f", t, T);
plot(T,sol(:,2));
axis([0 8 0 65])