% This function is a front-end or driver for the coupled oscillators problem.
function no_result = run_oscillators()
   
% Physical setup:  (shared using global variables)
global m1 m2 k1 k2 k3 L1 L2 L3
m1 = 1.0
m2 = 1.0
k1 = 1.0
k2 = 1.0
k3 = 0.3
L1 = 1.0 % Equilibrium value of x1 for spring 1.
L2 = 2.0 % Equilibrium value of x2 for spring 2.
L3 = 1.0 % Equilibrium value of x2-x1 for spring 3.

% Initial conditions:
   x1_initial = 1.2 * L1
   x2_initial = L2
   x1dot_initial = 0.0
   x2dot_initial = 0.0
   y0 = [x1_initial, x2_initial, x1dot_initial, x2dot_initial];
   
% Timing:
   time_step = 0.1
   total_time = 100
   N_steps = total_time/time_step
   
   % Create tabulated time values
   time = 0.0:time_step:total_time;
   
   % Solve the differential equation for y(time)
   y = rk2(y0, 0.0, time_step, N_steps, 'oscillators');

% Presentation of results:
   x1 = y(:,1);
   x2 = y(:,2);
   v1 = y(:,3);
   v2 = y(:,4);

   hold off
   plot(time,x1,'r')
   hold on
   plot(time,x2,'b')
   hold off
