function [t,y] = SMDRF(k, m, w, e, ti, tf)     % ‘SMRF’ for ‘Spring-Mass Rotating Frame’

% [t,y] = SDRF(104000,1/16,1,1,0,1000)

% Physical Parameters

M = 200;
C = 10;

% Initial Conditions

initialRad = 0.0;
initialRadVel = 0;
initialAngle = 0;
initialAngleVel = 0;

% Tolerances

tol = odeset('RelTol', 1e-8, 'AbsTol', 1e-9);

%     initial X, initial X', initial Y, initial Y'
ic = [initialRad; initialRadVel; initialAngle; initialAngleVel];
tspan = [ti tf];

% Equations of motion
% y(1) = X1
% y(2) = X2 = X1'
% y(3) = Y1
% y(4) = Y2 = Y1'

[t,y] = ode45(@odefun, tspan, ic, tol);

function dydt = odefun(t,y)
        dydt = zeros(4,1);
        dydt(1) = y(2);
        % dydt(2) = -(k/M)*y(1)+(m/M)*2*y(4)*w+(m/M)*(y(1)+e)*w^2; % No Damping
        dydt(2) = -(k/M)*y(1)-(C/M)*(y(2)-y(3)*w)+(m/M)*2*y(4)*w+(m/M)*(y(1)+e)*w^2; % Damping
        dydt(3) = y(4);
        % dydt(4) = -(k/M)*y(3)-(m/M)*2*y(2)*w+(m/M)*y(3)*w^2; %No damping
        dydt(4) = -(k/M)*y(3)-(C/M)*(y(1)+y(4)*w)-(m/M)*2*y(2)*w+(m/M)*y(3)*w^2; % Damping
end

%hold all
%plot(t,y(:,1),'-or')
%plot(t,y(:,3),'-ob')
%plot(t,(y(:,1).*y(:,1)+y(:,3).*y(:,3)).^0.5,'-og')

end