% Neil Edelman (110121860)

% plot the actual data
x = [.5:.5:7];
y = [5.4,-5.0,4.6,-4.4,3.9,-3.75,3.2,-3.25,2.75,-2.7,2.15,-2.05,1.95,-1.90];
figure(1);
plot(x, y, ' xr');
xlabel 'Oscillations, t/\tau';
ylabel 'Amplitude, A (m)';
title 'Unforced, Damped Harmonic Oscillator, Neil Edelman (110121860)'

% settings
a = 0;
b = 7;
res = 256;

dn = (b - a) / res;
% preallocate the solution matrix and a variable matrix
clear f;
clear x;
f(res + 1) = 0;
x = [a:dn:b];

% fill it in
f(1) = -6.17;
diff = 0;
n = a;
gamma = 0.42;
beta = 0;
omega2 = 1.141 .^ 2;
for i = 1:res
    k1 = -gamma * diff            - beta * diff .^ 2            - omega2 * f(i);
    k2 = -gamma * (diff + k1 / 2) - beta * (diff + k1 / 2) .^ 2 - omega2 * f(i);
    k3 = -gamma * (diff + k2 / 2) - beta * (diff + k2 / 2) .^ 2 - omega2 * f(i);
    k4 = -gamma * (diff + k3)     - beta * (diff + k3) .^ 2     - omega2 * f(i);
    diff = diff + dn * (k1 + 2 * k2 + 2 * k3 + k4) / 6;
    f(i + 1) = f(i) + diff;
    n = n + dn;
end
hold on;
plot(x, f, '- b');
hold off;
