dt=1e-4; %10^-3
m=1;
g=9.81;

A=m; %inertia, resistance to a change in velocity
B=4; %resistive force that pushes against motion
C=1; %restoring force, like a spring pushing back from being moved from rest
F=@(t) 0*sin(1.5*t); %external force

%calculates the next height/velocity/acceleration
nh=@(h,v) h+v*dt;
nv=@(v,a) v+a*dt;
na=@(h,v,t) (1/A)*(F(t)-B*v-C*h);

%initial conditions
h0=1;
v0=0;
a0=0;
Nt=500000;

%coding shit to prepare our 'table' of values
t=[0:dt:dt*Nt];
h=zeros(Nt+1,1);
v=zeros(Nt+1,1);
a=zeros(Nt+1,1);

%putting the initial conditions into the table
h(1)=h0;
v(1)=v0;
a(1)=a0;

%calculating the next height/velocity/acceleration at each time
for time=2:Nt+1
    h(time)=nh(h(time-1),v(time-1));
    v(time)=nv(v(time-1),a(time-1));
    a(time)=na(h(time),v(time),t(time));
end

%plot the result!
figure(1)
clf
plot(t,h,'k')
%the real solution we calced to compare (for A=m, B=0, C=0, F=-mg)
%hreal=(-g/2)*(t.^2)+v0*t+h0;
%hold on
%plot(t,hreal,'r')