// System component values
mass_1 = 1;
mass_2 =1;
spr_1 = 1;
spr_2 = 1;
spr_3 = 1;
d_1 = 1;
d_2 = 1;
d_3 = 1;
force = 1;

x1_0 = 0;		// initial conditions
v1_0 = 0;
x2_0 = 0;
v2_0 = 0;
X=[x1_0, v1_0, x2_0, v2_0]; //State Vector

// define the state matrix function
// the values returned are [x, v]
function foo=f(state,t)
	foo = [ state($, 2), ((-d_1-d_2)/mass_1)*state($,2)+((-spr_1-spr_2)/mass_1)*state($,1)+(d_2/mass_1)*state($,4)+(spr_2/mass_1)*state($,3), state($,4), ((-d_2-d_3)/mass_2)*state($,4)+((-spr_2-spr_3)/mass_2)*state($,3)+(d_2/mass_2)*state($,2)+(spr_2/mass_2)*state($,1)+ force/mass_2]; 
endfunction


// Set the time length and step size for the integration
steps = 100;
t_start = 3;
t_end = 20;
h = (t_end - t_start) / steps;


//
// Loop for integration
//
for i=1:steps,
	X = [X ; X($,:) + h*f(X, i*h)];
end
printf("The value at the end of first order integration is (x1, v1, x2, v2) =\n (%f, %f, %f, %f)\n", ...
	X($,1),	...
	X($,2),     ...
	X($,3),     ...
	X($,4));