//Scilab function for solving state equations
//By: Clayton Campagna, Sept, 16, 2004

// Get System component Values
mprintf("Enter the value of m1-->");	
[m1]=scanf("%f");
mprintf("Enter the value of m2");
[m2]=scanf("%f");
mprintf("Enter the value of Kd1");
[kd1]=scanf("%f");
mprintf("Enter the value of ks1");
[ks1]=scanf("%f");
mprintf("Enter the value of kd2");
[kd2]=scanf("%f");
mprintf("Enter the value of ks2");
[ks2]=scanf("%f");
mprintf("Enter the value of kd3");
[kd3]=scanf("%f");
mprintf("Enter the value of ks3");
[ks3]=scanf("%f");
mprintf("Enter the value of force");
[force]=scanf("%f");
mprintf("Enter the lower integration limit");
[t_start]=scanf("%f");
mprintf("Enter the upper integration limit");
[t_end]=scanf("%f");
mprintf("Enter the number of steps");
[steps]=scanf("%f");

x1=0;	//initial conditions
x2=0;
v1=0;
v2=0;
X=[x1, x2, v1, v2];

//define state matrix funcion
//the values returned are [x1,x2,v1,v2]

function foo=f(state, t)
	foo = [ state($, 3), 1*state($, 4), ((-ks1-ks2)/m1)*state($, 1) + (ks2/m1)*state($, 2) + ((-kd1-kd2)/m1)*state($, 3)	+ (kd2/m1)*state($, 4), (ks2/m2)*state($, 1) + ((-ks2-ks3)/m2)*	state($, 2) + (kd2/m2)*state($, 3) + ((-kd2-kd3)/m2)*state($, 4) + (force/m2)];
endfunction

//set the time length and step size for integration
h = (t_end - t_start) / steps;
t = [t_start];

//loop for integration
for i = 1:steps, 
	t = [t ; t($,:) + h];
	F1 = h * f(X($,:), t($,:));
	F2 = h * f(X($,:) + F1/2.0, t($,:) + h/2.0);
	F3 = h * f(X($,:) + F2/2.0, t($,:) + h/2.0);
	F4 = h * f(X($,:) + F3, t($,:) + h);
	X = [X ; X($,:) + (F1 + 2.0*F2 + 2.0*F3 + F4)/6.0];
end
X

//print some results to compare
mprintf("The value at the end of the first order integration is (x1, x2, v1, v2) = (%f, %f, %f, %f)\n", ...
	X($,1),	X($,2),	X($,3),	X($,4));