here's the code:
theta0=40;
v0=25;
deltat=0.01;
[x,y,t]=projectile_simple(v0, theta0, deltat);
%%finite diference
%forward
Vx_f=zeros(1,length(x));
Vy_f=zeros(1,length(x));
for jj=1:length(x)-1
Vx_f(jj)=(x(jj+1)-x(jj))/deltat;
Vy_f(jj)=(y(jj+1)-y(jj))/deltat;
end
%need correction for last index
Vx_f(end)=Vx_f(end-1);
Vy_f(end)=Vy_f(end-1);
figure
hold on
plot(t,Vy_f)
%backward
Vx_b=zeros(1,length(x));
Vy_b=zeros(1,length(x));
for k=2:length(x)
Vx_b(k)=(x(k-1)-x(k))/-deltat;
Vy_b(k)=(y(k-1)-y(k))/-deltat;
end
%need correction for last index
Vx_b(1)=Vx_b(2);
Vy_b(1)=Vy_b(2);
%deltat=x(2)-x(1);
plot(t,Vy_b)
%
%central
Vx_c=zeros(1,length(x));
Vy_c=zeros(1,length(x));
for m=2:length(x)-1
Vx_c(m)=((x(m+1)-x(m-1)))/(2*deltat);
Vy_c(m)=((y(m+1)-x(m-1)))/(2*deltat);
end
%need correction for last index
Vx_c(1)=Vx_c(2);
Vy_c(1)=Vy_c(2);
Vx_c(end)=Vx_c(end-1);
Vy_c(end)=Vy_c(end-1);
deltat=x(2)-x(1);
plot(t,Vy_c)
