x1_initial = 0;
x2_initial = 0;
v1_initial = 0;
v2_initial = 0;

Kd1 = 20;
Kd2 = 30;
Kd3 = 40;

Ks1 = 20;
Ks2 = 30;
Ks3 = 40;

F = 5;

M1 = 50;
M2 = 50;

start_time = 0;
end_time = 25;
steps = 1000;

h = (end_time - start_time) / steps;

t = [start_time];

A = [0, 0, -((Ks1 + Ks2)/M1), (Ks2/M2);
     0, 0, (Ks2/M1), -((Ks2 + Ks3)/M2);
     1, 0, -((Kd1 + Kd2)/M1), (Kd2/M2);
     0, 1, (Kd2/M1), -((Kd2 + Kd3)/M2)];
     
Bu = [0, 0, 0, (F/M2)];

X = [x1_initial, x2_initial, v1_initial, v2_initial];
X_new = [x1_initial, x2_initial, v1_initial, v2_initial];
printf("\n%f", t);

for i = 1:steps,

    F1 = h*(X*A + Bu);
    F2 = h*((X + F1/2)*A + Bu);
    F3 = h*((X + F2/2)*A + Bu);
    F4 = h*((X + F3)*A + Bu);
    X = [X + (F1 + 2*F2 + 2*F3 + F4)/6];
    X_new = [X_new ; X];
    t = t + h;
    printf("\n%f", t);
    
end

for j = 1:4,

    for i = 1:steps,

            printf("\n%f", X_new(i, j));

    end
    
end