k2 = 4.66;
k1 = k2;
m1 = 0.0917;
m2 = 0.0765;
number_of_ = 10;
compare = ones(40, 1000001);
n = 100000:100000:1000000;


A = [0, 1, 0, 0; ((-k1-k2)/m1), 0, (k2/m1), 0; 0, 0, 0, 1; (k2/m2), 0, (-k2/m2), 0];
for j = 1:length(n)
    y = zeros(4, n(j)+1);
    y(:, 1) = [100, 0, 50, 0];
    y_dot = zeros(4, n(j)+1);
    delta_t = 10/n(j);
    xi = 0:(10/n(j)):10;

%figure out how to assign these values to a matrix dependent on N
    for i = 2:length(xi)
        y(:, i) = y(:, i-1) + delta_t.*(A*y(:, i-1));
%        y_dot(:, i-1) = A*y(:, i-1);
    end
    compare(1+4*(j-1):4+4*(j-1), 1:length(y)) = y(:,:);

end


subplot(2,1,1);
plot(0:(10/n(1)):10, compare(1, 1:n(1)+1), '--y', 'Linewidth', 1.25);
hold on
plot(0:(10/n(2)):10, compare(5, 1:n(2)+1), '--m', 'Linewidth', 1.25);
hold on
plot(0:(10/n(3)):10, compare(9, 1:n(3)+1), '--c', 'Linewidth', 1.25);
hold on
plot(0:(10/n(4)):10, compare(13, 1:n(4)+1), '--r', 'Linewidth', 1.25);
hold on
plot(0:(10/n(5)):10, compare(17, 1:n(5)+1), '--g', 'Linewidth', 1.25);
hold on
plot(0:(10/n(6)):10, compare(21, 1:n(6)+1), '--b', 'Linewidth', 1.25);
hold on
plot(0:(10/n(7)):10, compare(25, 1:n(7)+1), '--y', 'Linewidth', 1.25);
hold on
plot(0:(10/n(8)):10, compare(29, 1:n(8)+1), '--m', 'Linewidth', 1.25);
hold on
plot(0:(10/n(9)):10, compare(33, 1:n(9)+1), '--c', 'Linewidth', 1.25);
hold on
plot(0:(10/n(10)):10, compare(37, 1:n(10)+1), '--r', 'Linewidth', 1.25);
hold on
legend('N = 100000', 'N = 200000','N = 300000', 'N = 400000', 'N = 500000', 'N = 600000', 'N = 700000', 'N = 800000', 'N = 900000', 'N = 1000000');
title('Approximating position')
xlabel('Time');
ylabel('Position of x1(t)');
xlim([0 10]);
grid on;


subplot(2,1,2);
plot(0:(10/n(1)):10, compare(3, 1:n(1)+1), '--y', 'Linewidth', 1.25);
hold on;
plot(0:(10/n(2)):10, compare(7, 1:n(2)+1), '--m', 'Linewidth', 1.25);
hold on;
plot(0:(10/n(3)):10, compare(11, 1:n(3)+1), '--c', 'Linewidth', 1.25);
hold on;
plot(0:(10/n(4)):10, compare(15, 1:n(4)+1), '--r', 'Linewidth', 1.25);
hold on;
plot(0:(10/n(5)):10, compare(19, 1:n(5)+1), '--g', 'Linewidth', 1.25);
hold on;
plot(0:(10/n(6)):10, compare(23, 1:n(6)+1), '--b', 'Linewidth', 1.25);
hold on;
plot(0:(10/n(7)):10, compare(27, 1:n(7)+1), '--y', 'Linewidth', 1.25);
hold on;
plot(0:(10/n(8)):10, compare(31, 1:n(8)+1), '--m', 'Linewidth', 1.25);
hold on;
plot(0:(10/n(9)):10, compare(35, 1:n(9)+1), '--c', 'Linewidth', 1.25);
hold on;
plot(0:(10/n(10)):10, compare(39, 1:n(10)+1), '--r', 'Linewidth', 1.25);
hold on;
legend('N = 100000', 'N = 200000','N = 300000', 'N = 400000', 'N = 500000', 'N = 600000', 'N = 700000', 'N = 800000', 'N = 900000', 'N = 1000000');
title('Approximating position')
xlabel('Time');
ylabel('Position of x2(t)');
xlim([0 10]);
grid on;

eigen = eig(A)