Exam3_3

function exam3_3()
    % Initialize step-size variables
    h = 0.1;
    omega_plot = (0:0.1:7)';
    t = (0:0.1:7)';
    N = length(t);
    % each amplitudes
    amplitude1 = zeros(N-1,1);
    amplitude2 = zeros(N-1,1);
    amplitude3 = zeros(N-1,1);
    amplitude4 = zeros(N-1,1);
    amplitude5 = zeros(N-1,1);
    % max amplitude
    max_amplitude = zeros(N-1,1);
    % Initialize derivative function
    dz1 = @(t,x1,x2,omega)(-20*x1 + 10*x2 - ((omega^2)*cos(omega*t)));
    dz2 = @(t,x1,x2,x3,omega)(10*x1 - 20*x2 + 10*x3 - ((omega^2)*cos(omega*t)));
    dz3 = @(t,x2,x3,x4,omega)(10*x2 - 20*x3 + 10*x4 - ((omega^2)*cos(omega*t)));
    dz4 = @(t,x3,x4,x5,omega)(10*x3 - 20*x4 + 10*x5 - ((omega^2)*cos(omega*t)));
    dz5 = @(t,x4,x5,omega)(10*x4 - 20*x5 - ((omega^2)*cos(omega*t)));
    dx1 = @(t,z1) z1;
    dx2 = @(t,z2) z2;
    dx3 = @(t,z3) z3;
    dx4 = @(t,z4) z4;
    dx5 = @(t,z5) z5;
    % Iterate, computing each K value in turn, then the i+1 step values
    counter = 1;
    for omega = 0:0.1:7
        % Initialize vectors
        x1 = zeros(N,1); x2 = zeros(N,1); x3 = zeros(N,1); x4 = zeros(N,1); x5 = zeros(N,1);
        z1 = zeros(N,1); z2 = zeros(N,1); z3 = zeros(N,1); z4 = zeros(N,1); z5 = zeros(N,1);
        % Initialize K vectors
        kx1 = zeros(1,4)'; kx2 = zeros(1,4)'; kx3 = zeros(1,4)'; kx4 = zeros(1,4)'; kx5 = zeros(1,4)';
        kz1 = zeros(1,4)'; kz2 = zeros(1,4)'; kz3 = zeros(1,4)'; kz4 = zeros(1,4)'; kz5 = zeros(1,4)';
        b = [1 2 2 1]'; % RK4 coefficients
        for i = 1:(N-1) %% for t
            kz1(1) = dz1(t(i),x1(i),x2(i),omega);
            kz2(1) = dz2(t(i),x1(i),x2(i),x3(i),omega);
            kz3(1) = dz3(t(i),x2(i),x3(i),x4(i),omega);
            kz4(1) = dz4(t(i),x3(i),x4(i),x5(i),omega);
            kz5(1) = dz5(t(i),x4(i),x5(i),omega);
            kx1(1) = z1(i); kx2(1) = z2(i); kx3(1) = z3(i); kx4(1) = z4(i); kx5(1) = z5(i);

            kz1(2) = dz1(t(i) + (h/2),x1(i)+(h/2)*kx1(1),x2(i)+(h/2)*kx2(1),omega);
            kz2(2) = dz2(t(i) + (h/2),x1(i)+(h/2)*kx1(1),x2(i)+(h/2)*kx2(1),x3(i)+(h/2)*kx3(1),omega);
            kz3(2) = dz3(t(i) + (h/2),x2(i)+(h/2)*kx2(1),x3(i)+(h/2)*kx3(1),x4(i)+(h/2)*kx4(1),omega);
            kz4(2) = dz4(t(i) + (h/2),x3(i)+(h/2)*kx3(1),x4(i)+(h/2)*kx4(1),x5(i)+(h/2)*kx5(1),omega);
            kz5(2) = dz5(t(i) + (h/2),x4(i)+(h/2)*kx4(1),x5(i)+(h/2)*kx5(1),omega);
            kx1(2) = z1(i) + (h/2)*kz1(1);
            kx2(2) = z2(i) + (h/2)*kz2(1);
            kx3(2) = z3(i) + (h/2)*kz3(1);
            kx4(2) = z4(i) + (h/2)*kz4(1);
            kx5(2) = z5(i) + (h/2)*kz5(1);

            kz1(3) = dz1(t(i) + (h/2),x1(i)+(h/2)*kx1(2),x2(i)+(h/2)*kx2(2),omega);
            kz2(3) = dz2(t(i) + (h/2),x1(i)+(h/2)*kx1(2),x2(i)+(h/2)*kx2(2),x3(i)+(h/2)*kx3(2),omega);
            kz3(3) = dz3(t(i) + (h/2),x2(i)+(h/2)*kx2(2),x3(i)+(h/2)*kx3(2),x4(i)+(h/2)*kx4(2),omega);
            kz4(3) = dz4(t(i) + (h/2),x3(i)+(h/2)*kx3(2),x4(i)+(h/2)*kx4(2),x5(i)+(h/2)*kx5(2),omega);
            kz5(3) = dz5(t(i) + (h/2),x4(i)+(h/2)*kx4(2),x5(i)+(h/2)*kx5(2),omega);
            kx1(3) = z1(i) + (h/2)*kz1(2);
            kx2(3) = z2(i) + (h/2)*kz2(2);
            kx3(3) = z3(i) + (h/2)*kz3(2);
            kx4(3) = z4(i) + (h/2)*kz4(2);
            kx5(3) = z5(i) + (h/2)*kz5(2);

            kz1(4) = dz1(t(i) + h,x1(i)+h*kx1(3),x2(i)+h*kx2(3),omega);
            kz2(4) = dz2(t(i) + h,x1(i)+h*kx1(3),x2(i)+h*kx2(3),x3(i)+h*kx3(3),omega);
            kz3(4) = dz3(t(i) + h,x2(i)+h*kx2(3),x3(i)+h*kx3(3),x4(i)+h*kx4(3),omega);
            kz4(4) = dz4(t(i) + h,x3(i)+h*kx3(3),x4(i)+h*kx4(3),x5(i)+h*kx5(3),omega);
            kz5(4) = dz5(t(i) + h,x4(i)+h*kx4(3),x5(i)+h*kx5(3),omega);
            kx1(4) = z1(i) + h*kz1(3);
            kx2(4) = z2(i) + h*kz2(3);
            kx3(4) = z3(i) + h*kz3(3);
            kx4(4) = z4(i) + h*kz4(3);
            kx5(4) = z5(i) + h*kz5(3);

            z1(i+1) = z1(i) + (h/6)*sum(b.*kz1);
            z2(i+1) = z2(i) + (h/6)*sum(b.*kz2);
            z3(i+1) = z3(i) + (h/6)*sum(b.*kz3);
            z4(i+1) = z4(i) + (h/6)*sum(b.*kz4);
            z5(i+1) = z5(i) + (h/6)*sum(b.*kz5);
            x1(i+1) = x1(i) + (h/6)*sum(b.*kx1);
            x2(i+1) = x2(i) + (h/6)*sum(b.*kx2);
            x3(i+1) = x3(i) + (h/6)*sum(b.*kx3);
            x4(i+1) = x4(i) + (h/6)*sum(b.*kx4);
            x5(i+1) = x5(i) + (h/6)*sum(b.*kx5);
        end
        amplitude1(counter) = max(x1) - min(x1);
        amplitude2(counter) = max(x2) - min(x2);
        amplitude3(counter) = max(x3) - min(x3);
        amplitude4(counter) = max(x4) - min(x4);
        amplitude5(counter) = max(x5) - min(x5);
        temp = [amplitude1(counter); amplitude2(counter); amplitude3(counter); amplitude4(counter); amplitude5(counter);];
        max_amplitude(counter) = max(temp);
        counter = counter + 1;
    end
    plot(omega_plot,max_amplitude,'b-');
end