%solve dv/dt
n=1;

function xdot = f(x,t)
  c1 = 0.25;
  m = 0.5;
  % x(2)=x,  x(1)=v
  xdot(1) = x(2);
  xdot(2) = -(c1/m)*x(2);
endfunction

function ydot = g(y,t)
  c2 = 0.50;
  m = 0.5;
  % y(2)=y, y(1)=v
  ydot(1) = y(2);
  ydot(2) = -(c2/m)*y(2);
endfunction

function zdot = h(z,t)
  g = 9.81;
  c3 = 0.9;
  m = 0.5;
  %z(1) = z, z(2) = zdot(1) = vz, zdot(2) = az
  zdot(1) = z(2);
  zdot(2) = ((-c3*z(2))-m*g)/m;
endfunction

g=9.81;
m=0.5;
c3=0.9;
c2=0.5;
c1=0.25;
k=m/c3;
vx0 = 25; vy0 = 15; vz0 = 10;
x0 = 0; y0 = 0; z0 = 0;


x0vx0=[0,25];
y0vy0=[0,15];
z0vz0=[,10];
t=linspace(0,2.5,50);


xt=vx0*e.^(-c1*t/m);
yt=vy0*e.^(-c2*t/m);
xt=vx0*e.^(-t/k)-k*g*e.^(-t/k)-k*g;;

xxt=lsode("f",x0vx0,t);
yyt=lsode("g",y0vy0,t);
zzt=lsode("h",z0vz0,t);


figure(1);
plot3(xxt(:,1),yyt(:,1),zzt(:,1),'-r','linewidth',2);
grid on;
title(" Projectile Motion ", 'fontsize',14);
xlabel("x position (m)",'fontsize',14);
ylabel("y position (m)",'fontsize',14);
zlabel("z position (m)",'fontsize',14);
axis([-1 5 5 1 0 20]);



figure(2);
plot3(xt(:,1),yt(:,1),zt(:,1),'-r','linewidth',2);
grid on;
title(" Projectile Motion ", 'fontsize',14);
xlabel("x position (m)",'fontsize',14);
ylabel("y position (m)",'fontsize',14);
zlabel("z position (m)",'fontsize',14);
axis([-1 5 5 1 0 20]);