%Demo #1 Dynamics 2030, Spring 2019
n= 1000; h = 2*pi/n;     % h = delta t
k = 1;   m = 1;          % parameters
x = 1;   v = 0;          % initial conditions
xarray = x;  varray = v; % Start at the beginning

for i = 1:n
  xdot   = v;            % You could put most any  
  vdot   = -k*x/m;       % problem from Math 2930 here.
  x      = x + xdot * h; % THESE TWO LINES ARE THE...
  v      = v + vdot * h; % HEART AND SOLE OF THIS CODE.
  xarray = [xarray;x]; 
  varray = [varray;v];
end

plot (xarray, varray); axis('equal'); grid on;shg
disp(['Numerical error = ', num2str(x(end)-1)]);  
