param.m = 1;                            % Mass [kg]
param.k = 1;                            % Spring rate [N/m]
param.wn = sqrt(param.k/param.m);       % (Circular) natural frequency [rad/s]

x0 = 1;              % Position [m]
x0_dot = 0;          % Velocity [m/s]
z0 = [x0 ; x0_dot];  % Initial condition on state vector

tend = 4*round((2*pi)/param.wn);    % Simulate for approximately 4 periods [s]
dt = round((2*pi)/param.wn)/1000;   % Time step [s]
tsim = 0:dt:tend;                   % Vector of simulation times [s]

param.zeta = 0;     % Damping ratio

[t,z] = ode45(@(t,z) odesys( t,z,param ),tsim,z0);
x(:,1) = z(:,1);    % Save solution for future plotting.

plot(t,x(:,1),'linewidth',2);grid on;hold all

    ylim([-1 1])
    xlabel('Time [s]');
    ylabel('Position [m]')
    title('Free Response: Undamped 2nd Order System');

zetavals = [0 1/sqrt(2) 1 1.5]; % Damping ratios to consider.

case1 = 'Undamped: zeta = 0' ;
case2 = 'Underdamped: zeta = 0.707' ;
case3 = 'Critically Damped: zeta = 1' ;
case4 = 'Overdamped: zeta = 1.5' ;

for i = 2: length(zetavals)
    
    param.zeta = zetavals(i); % Set current damping ratio.
    [t,z] = ode45(@(t,z) odesys( t,z,param ),tsim,z0); % Compute free response.
    x(:,i) = z(:,1); % Save solution for future plotting.

end;

figure      % Plot results

    plot(t,x,'linewidth',2);grid on
    xlabel('Time [s]');
    ylabel('Position [m]')
    title('Free Response of 2nd Order System for Various Damping Ratios');
    legend(case1,case2,case3,case4,0)
    
 
zeta_range = 0:0.01:1.5;          % Damping ratio: undamped to overdamped.
S = zeros(2,length(zeta_range));  % Initialize characteristic roots array.

for i = 1:length(zeta_range)
    
    S(:,i) = eig([0 1; -param.wn^2 -2*zeta_range(i)*param.wn]); % Characteristic roots.

end;

figure      % Plot characteristic roots in complex plane.

    plot(S,'bx','linewidth',2,'markersize',5); grid on;
    set(gca,'YAxisLocation','right')
    xlabel('Real Axis','Position',[-1.5 -0.1 1])
    ylabel('Imaginary Axis')
    title('Loci of Characteristic Roots in Complex Plane for Various Damping Ratios')

type odesys.m