// System component values
m1 = 100;
m2 = 38;
F = 1000;
kd1 = 1.5;
kd2 = 5.1;
kd3 = 0.5;
ks1 = 3.2;
ks2 = 1.8;
ks3 = 7.2;

x1 = 0;		// initial conditions
x2 = 0;
v1 = 0;
v2 = 0;

X=[x1, v1, x2, v2];

// define the state matrix function
// the values returned are [x1, v1, x2, v2]
function foo=f(state,t)
	foo = [ state($, 2), state($, 2)*((-kd1-kd2)/m1)+state($, 1)*((-ks1-ks2)/m1)+...
		state($, 4)*(kd2/m1)+state($, 3)*(ks2/m1),state($, 4),state($, 4)*...
		((-kd2-kd3)/m2)+state($, 3)*((-ks2-ks3)/m2)+state($, 2)*(kd2/m2)+...
                state($, 1)*(ks2/m2)+(F/m2)];
endfunction


// Set the time length and step size for the integration
steps = 1000;
t_start = 0;
t_end = 4.5;
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) = (%f, %f, %f, %f)\n", ...
	X($,1),	...
	X($,2), ...
	X($,3), ...
	X($,4));