(**************************************)
(*Initialisation*)(* 1D motion with drag force and constant P *)
Deltat=0.01;tini=0.001;tlast=1.0;
P=1;m=1.0;y[0]=10;
g[i_]:=P/(2*m*y[i]);
(**************************************)

(**************************************)
(* Calculation *)
Ni=IntegerPart[(tlast-tini)/Deltat];
For[i=1,i≤Ni,i++,
y[i]=y[i-1]+g[i-1]*Deltat;
t[i]=tini+i*Deltat;
];(*end for *)
tty= Table[{t[i],y[i]},{i,0,Ni}];
(**************************************)

(* Output *)
(**************************************)
ListPlot[tty,PlotRange->All]
(**************************************)

(**************************************)
(*Initialisation*)
Deltat=0.01;tini=0.001;tlast=50.0;
P=1;m=1.0;y[0]=10;
c=1;A=1;rho=1;

g[i_]:=P/(2*m*y[i]) - (rho*c*A*y[i]/(2*m));
(**************************************)

(**************************************)
(* Calculation *)
Ni=IntegerPart[(tlast-tini)/Deltat];
For[i=1,i≤Ni,i++,
y[i]=y[i-1]+g[i-1]*Deltat;
t[i]=tini+i*Deltat;
];(*end for *)
tty= Table[{t[i],y[i]},{i,0,Ni}];
(**************************************)

(* Output *)
(**************************************)
ListPlot[tty,PlotRange->All]
(**************************************)

(**************************************)
(*Initialisation*)
Deltat=0.1;tini=0.0;tlast=10.0;
g=9.81;
vy[0]=0;y[0]=0;t[0]=0;y[1]=y[0]+Deltat*vy[0];
F[i_]:=g;
(**************************************)

(**************************************)
(* Calculation *)
Ni=IntegerPart[(tlast-tini)/Deltat];

For[i=1,i≤Ni,i++,
t[i]=tini+i*Deltat;
vy[i]=vy[i-1]+F[i-1]*Deltat;
y[i+1] = y[i] + vy[i]*Deltat;
];(*end for *)
ttvy= Table[{t[i],vy[i]},{i,0,Ni}];
tty= Table[{t[i],y[i]},{i,0,Ni}];
(**************************************)

(* Output *)
(**************************************)
ListPlot[ttvy,PlotRange->All,AxesLabel->{"t","vy "}]
ltty=ListPlot[tty,PlotRange->All,AxesLabel->{"t","y"}];
ptty=Plot[0.5*g*t^2,{t,tini,tlast},PlotStyle->Red];
Show[ltty,ptty]

(**************************************)

(**************************************)
(*Initialisation*)(* 2D Projectile without drag force  *)
Deltat=0.1;tini=0.0;tlast=10.0;
g=9.81;
vy[0]=20.0;y[0]=0;t[0]=0;y[1]=y[0]+Deltat*vy[0];
vx[0]=1;x[0]=0;x[1]=x[0]+Deltat*vx[0];

F[i_]:=-g;  Fx[i_]:=0;
(**************************************)

(**************************************)
(* Calculation *)
Ni=IntegerPart[(tlast-tini)/Deltat];

For[i=1,i≤Ni,i++,
t[i]=tini+i*Deltat;
vy[i]=vy[i-1]+F[i-1]*Deltat;
y[i+1] = y[i] + vy[i]*Deltat;

vx[i]=vx[i-1]+Fx[i-1]*Deltat;
x[i+1] = x[i] + vx[i]*Deltat;

];(*end for *)
ttvy= Table[{t[i],vy[i]},{i,0,Ni}];
tty= Table[{t[i],y[i]},{i,0,Ni}];

ttvx= Table[{t[i],vx[i]},{i,0,Ni}];
ttx= Table[{t[i],x[i]},{i,0,Ni}];
txy=Table[{x[i],y[i]},{i,0,Ni}];
(**************************************)

(* Output *)
(**************************************)
ListPlot[ttvy,PlotRange->All,AxesLabel->{"t","vy "}];
ltty=ListPlot[tty,PlotRange->All,AxesLabel->{"t","y"}];
ptty=Plot[0.5*g*t^2,{t,tini,tlast},PlotStyle->Red];

lttx=ListPlot[ttx,PlotRange->All,AxesLabel->{"t","x"}];

ltxywod=ListPlot[txy,PlotRange->All,AxesLabel->{"x","y"}];
Show[ltxywod]
(**************************************)

txy

Manipulate[
ListPlot[{{x[i],y[i]}},PlotRange->{{0,x[Ni]},{0,y[Ni]}},PlotStyle->{PointSize->0.04}],{i,0,Ni,1}]

txy

lasttxy=Last[txy]

xlast=lasttxy[[1]]

ylast=lasttxy[[2]]

Manipulate[
ListPlot[{{x[i],y[i]}},PlotRange->{{0,xlast},{0,ylast}},PlotStyle->{PointSize->0.04}],{i,0,Ni,1}]

exam={{a11,a12},{a21,a22},{a31,a32}}

lastexam=Last[exam]

lastexam[[1]]

lastexam[[2]]

t=1;Plot[Sin[3x-0.1*t],{x,-Pi,Pi},PlotLabel->{"t=",t}];
t=2;Plot[Sin[3x-0.1*t],{x,-Pi,Pi},PlotLabel->{"t=",t}];
t=3;Plot[Sin[3x-0.1*t],{x,-Pi,Pi},PlotLabel->{"t=",t}];

Manipulate[
Plot[Sin[3x-0.1*t],{x,-Pi,Pi},PlotLabel->{"t=",t}]
,{t,0,10,1}]

(* the following code is for a 2D projectile motion with drag force  *)

(**************************************)
(*Initialisation*)
Deltat=0.05;tini=0.0;tlast=70;
g=9.81;m=1.0;
rho=1.0;A=10.0^(-4);Cd=100.0;
B2=0.5*Cd*rho*A;

vy[0]=350;y[0]=0;t[0]=0;y[1]=y[0]+Deltat*vy[0];
vx[0]=1;x[0]=0;x[1]=x[0]+Deltat*vx[0];

F[i_]:=-g;  Fx[i_]:=0;
(**************************************)

(**************************************)
(* Calculation *)
Ni=IntegerPart[(tlast-tini)/Deltat];

For[i=1,i≤Ni,i++,
t[i]=tini+i*Deltat;

v[i-1]=Sqrt[vx[i-1]^2+vy[i-1]^2];

vy[i]=vy[i-1]+F[i-1]*Deltat - B2*vy[i-1]*v[i-1]*Deltat/m;
y[i+1] = y[i] + vy[i]*Deltat;

vx[i]=vx[i-1]+Fx[i-1]*Deltat- B2*vx[i-1]*v[i-1]*Deltat/m;

x[i+1] = x[i] + vx[i]*Deltat;

];(*end for *)
ttvy= Table[{t[i],vy[i]},{i,0,Ni}];
tty= Table[{t[i],y[i]},{i,0,Ni}];

ttvx= Table[{t[i],vx[i]},{i,0,Ni}];
ttx= Table[{t[i],x[i]},{i,0,Ni}];
txy=Table[{x[i],y[i]},{i,0,Ni}];
(**************************************)

(* Output *)
(**************************************)
ListPlot[ttvy,PlotRange->All,AxesLabel->{"t","vy "}];
ltty=ListPlot[tty,PlotRange->All,AxesLabel->{"t","y"}];
ptty=Plot[0.5*g*t^2,{t,tini,tlast},PlotStyle->Red];

lttx=ListPlot[ttx,PlotRange->All,AxesLabel->{"t","x"}];

ltxy=ListPlot[txy,PlotRange->{{0,10},{-500,500}},AxesLabel->{"x","y"},PlotStyle->Red];
Show[ltxy(*,ltxywod*)]
(**************************************)

(* The following code is for a 2D projectile motion with (i) drag force, (ii) height-dependence of density of air  *)

(**************************************)
(*Initialisation*)
Deltat=0.05;tini=0.0;tlast=100;
g=9.81;m=1.0;
rho=1.0;A=10.0^(-4);Cd=100.0;
B2=0.5*Cd*rho*A;
y0=1000;
a=6.5*10^(-3);T0=300;alpha=2.5;

vy[0]=1000;y[0]=0;t[0]=0;y[1]=y[0]+Deltat*vy[0];
vx[0]=100;x[0]=0;x[1]=x[0]+Deltat*vx[0];

F[i_]:=-g;  Fx[i_]:=0;
(**************************************)

(**************************************)
(* Calculation *)
Ni=IntegerPart[(tlast-tini)/Deltat];

For[i=1,i≤Ni,i++,
t[i]=tini+i*Deltat;

v[i-1]=Sqrt[vx[i-1]^2+vy[i-1]^2];

imodel=2;
(*** choosing model of air density ************)
If[imodel==1,(*isothermal *) 
ratiorho[i-1]=Exp[-y[i-1]/y0]
];(*end if *)

If[imodel==2,(* adiabatic *) 
ratiorho[i-1]=(1 - a*y[i-1]/T0)^alpha;
];(*end if *)
(*** choosing model of air density ************)

vy[i]=vy[i-1]+F[i-1]*Deltat - (ratiorho[i-1])*B2*vy[i-1]*v[i-1]*Deltat/m;
y[i+1] = y[i] + vy[i]*Deltat;

vx[i]=vx[i-1]+Fx[i-1]*Deltat- (ratiorho[i-1])*B2*vx[i-1]*v[i-1]*Deltat/m;

x[i+1] = x[i] + vx[i]*Deltat;

];(*end for *)
ttvy= Table[{t[i],vy[i]},{i,0,Ni}];
tty= Table[{t[i],y[i]},{i,0,Ni}];

ttvx= Table[{t[i],vx[i]},{i,0,Ni}];
ttx= Table[{t[i],x[i]},{i,0,Ni}];
txy=Table[{x[i],y[i]},{i,0,Ni}];
(**************************************)

(* Output *)
(**************************************)
ListPlot[ttvy,PlotRange->All,AxesLabel->{"t","vy "}];
ltty=ListPlot[tty,PlotRange->All,AxesLabel->{"t","y"}];
ptty=Plot[0.5*g*t^2,{t,tini,tlast},PlotStyle->Red];

lttx=ListPlot[ttx,PlotRange->All,AxesLabel->{"t","x"}];

ltxy=ListPlot[txy,PlotRange->All,AxesLabel->{"x","y"},PlotStyle->Red];
Show[ltxy(*,ltxywod*)]
(**************************************)