//This proggy solves second order diff-eqs maybe even more.  Just be sure to change the variables.
//The below is for AP 5.9

ks1 = 100; 
ks2 = 100; 
kd = 10; 
m1 = 1;
m2 = 1;
F = 1

// initial conditions
x1 = 0; 
v1 = 0;
x2 = 0;
v2 = 0;
X = [x1, v1, x2, v2];	// Be sure to note the order of x1,v1,...etc



// define the state matrix function
function foo=f(state,t)
x1 = state($, 1);
v1 = state($, 2);
x2 = state($, 3);
v2 = state($, 4);

foo = [ v1, -v1*(kd/m1)+v2*(kd/m1)-x1*(ks1/m1)+x2*(ks1/m1)+(F/m1),v2,v1*(kd/m2)-v2*(kd/m2)+x1*(ks1/m1)-x2*(ks1/m2)-x2*(ks2/m2)];

endfunction


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

// Loop for integration
for i=1:steps,
	t = [t ; t_start + i*h];
	X = [X ; X($,:) + h*f(X, t_start + i*h)];
		time = i * h;
  		if time == 0 |time == 1 | time == 2 | time == 3 | time == 4 | time == 5...
 				| time == 6 | time == 7 | time == 8 | time == 9 | time == 10
 			printf("%f  %f\n",X($,1),X($,2));
 		end
end

plot2d(t,X,[9,3,6,5], leg="x1@v1@x2@v2");


//print_steps = 10;
//for i = 0:print_steps;
//	j = i * (steps / print_steps);
//	time = j * h + t_start;
//	printf("The value at t=%f s is (x, v) = (%f, %f)\n", time, X(j + 1,1), X(j +1,2));
//end