% This is a good idea if you are going to run this script many times
clear ;

% This is the system of units we are using

M = 1. ;
K = 1. ;
x0 = 1. ;
w0 = sqrt(K/M) ;

% We will take the initial velocity to be zero
% This is the initial velocity in units of  x0*w0
% i.e. this is vbar = v/ (x0*w0)
v0 = 0. ; 


% This defines the times step in units of w0
dt = 0.001 ;
% This defines the final time and initial time
tmin =0. ;
tmax =10.; 
% This defines an array t(1) = tmin , t(1) = tmin + dt, t(2) =tmin + 2dt,
% t(3) = tmin + 3*dt, ...
t = [tmin:dt:tmax]  ;

% We will store our solution in these arrays
% x(1) =0 , x(2) = 0, x(3) = 0...
% Later we will fill these arrays up with the actual solution
x = zeros(1, length(t) ) ; 
v = zeros(1, length(t) ) ;

% Later you might want to use this array to store the work done 
% by friction
Wfr = zeros(1, length(t) ) ;


xn = x0 ;
vn = v0 ;
tn = tmin ;

% Store this inital value of x and the initial value of velocity
x(1) = xn;
v(1) = vn;



for i = 2:length(t) 
   
    F = -K*xn ;
    
    % This step determines the n+1 values of xn, vn, tn
    xn=xn+ vn*dt ;
    vn=vn + F/M*dt ;
    tn=tn + dt;
    
           
    % Store the values of the  xn in the appropriate time slot.
    x(i) = xn ;
    v(i) = vn ;

    
end 

% This shows how to plot a function of the array t
% after this line
% xsol(1) = x0*cos(w0*t(1)), xsol(2)=x0*cos(w0*t(2)), ....
% This should only agree with the numerical results for B=0
xsol = x0*cos(w0*t) ;

% Uncomment this line to plot the  x versus t , 
% and the zero damping solution x= x0*cos(w0*t)
plot(t, x, t, xsol) ;