% This program also needs the function
% force.m , which must be in the same directory
% as this program
%
% Modify from here -->
N = 2;      % nr of balls
m = ...;    % kg
k = ...;    % N/m
q = 1.0;
d = ...;    % m
v0 = ...;   % m/s
time = ...; % s
dt = ...;   %s
% Base variables
n = ceil(time/dt);
x = zeros(n,N);
v = zeros(n,N);
F = zeros(n,N);
t = zeros(n,1);
% Initial conditions, equally spaced
for j = 1:N
    x(1,j) = d*(j-1);
end
v(1,1) = v0;
for i = 1:n-1
  % Find force on each block, j
  % First, force from block to the left
  for j = 2:N
    dx = x(i,j) - x(i,j-1);
    F(i,j) = F(i,j) + force(dx,d,k,q);
  end
  % Second, force from block to the right
  for j = 1:N-1
    dx = x(i,j+1) - x(i,j);
    F(i,j) = F(i,j) - force(dx,d,k,q);
  end
  % Euler-Cromer step
  for j = 1:N
    a = F(i,j)/m;
    v(i+1,j) = v(i,j) + a*dt;
    x(i+1,j) = x(i,j) + v(i+1,j)*dt;
  end
  %
  % The Euler-Cromer step above can also be vectorized
  % through the following implementation (which is faster)
  % Euler-Cromer vectorized step
  % a = F(i,:)/m;
  % v(i+1,:) = v(i,:) + a*dt;
  % x(i+1,:) = x(i,:) + v(i,:)*dt;
  %
  t(i+1) = t(i) + dt;
end
% Plot results
for j = 1:N
  plot(t,v(:,j));
  if j==1
    hold on;
  end
  if j==N
    hold off
  end
end
xlabel('t [s]');
ylabel('v [m/s]');
