%% Simulate a damped nonlinear spring, mass system
% Initiate the script
clc; clear; close all

%% Define the parameters
m = 1;      % kg
f = 0.6;    % Ns/m
k1 = 1e2;   % N/m
k3 = 5;     % N/m^3

% Initial conditions
IC1 = [1e-3; 0];
IC2 = [10; 0];

tt = linspace(0,3*pi,2e3);  % define the time vector - seconds
%% Simulate
% First IC
opts = odeset('RelTol',1e-4,'AbsTol',1e-6);
[t_out1,z_out1] = ode45(@(t,z)  Damped_nonlin_sys(t,z,m,f,k1,k3), tt, IC1,opts);
% Second IC
[t_out2,z_out2] = ode45(@(t,z)  Damped_nonlin_sys(t,z,m,f,k1,k3), tt, IC2,opts);

%% Data visualization  
% plot the normalized position vs. time
hf = figure(1); clf
set(hf,'units','Normalized','position',[0 0 1 1])
h = plot(tt,[z_out2(:,1)/max(z_out2(:,1)) z_out1(:,1)/max(z_out1(:,1))]);
    set(h(1),'linewidth',2);
    set(h(2),'linewidth',5);
    grid minor
    axis tight
    hl = legend(h,'IC_2','IC_1');
    xlabel('t(s)');
    ylabel('normalized x')
    set(gca,'FontSize',18)
    set(gcf,'color','w')
% plot the normalized position vs. normalized velocity    
hf = figure(2); clf
set(hf,'units','Normalized','position',[0 0 1 1])
subplot(1,2,1)
plot(z_out1(:,2)/max(z_out1(:,2)), z_out1(:,1)/max(z_out1(:,1)),'linewidth',3);
    axis equal
    grid minor
    title('IC_1');
    xlabel('normalized dx/dt');
    ylabel('normalized x')
    set(gca,'FontSize',18)
    set(gcf,'color','w')
subplot(1,2,2)
plot(z_out2(:,2)/max(z_out2(:,2)), z_out2(:,1)/max(z_out2(:,1)),'linewidth',3);
    axis equal
    grid minor
    title('IC_2');
    xlabel('normalized dx/dt');
    ylabel('normalized x')
    set(gca,'FontSize',18)
    
% %% Functions
% function dzdt = Damped_nonlin_sys(t,z,m,f,k1,k3)
% %Damped_nonlin_sys computes the derivative of the states z
% dzdt(1,1) = z(2);
% dzdt(2,1) = -(f*z(2)+k1*z(1)+k3*z(1)^3)/m;
% end