%% Problem 1
% 1 b)
syms l k l0
k(l) = (k/3)*((l-(l0^3/l^2))/(l-l0));
T = taylor(k(l), 'ExpansionPoint', l0,'Order', 3);

% 1c)
k0 = 16;
l0 = .2;

syms k(l)

k(l) = (k/3)*((l-(l0^3/l^2))/(l-l0));
l = linspace(0,0.5);
linear = k0*(l-l0);
nonlinear = -(k0/3)*((l-(l0^3./l.^2))./(l-l0));

figure(1)
subplot(211);
plot(linear,nonlinear,'LineWidth',1.5,'DisplayName',"Elastomer");
title('Problem 1: Question 1C');
xlabel('Distance');
ylabel('Elastomer Force');

subplot(212);
plot(l,linear,'LineWidth',1.5,'DisplayName',"Linear");
xlabel('Distance');
ylabel('Linear Force');

%2b)
l0=0.2;
m=1;
t0=0;
tf=6;

lt0 = 0.1;
v0=0;

tspan=[t0 tf];
u0 = [lt0 v0];
f1=@(t,u) [u(2), -1*((k0)*(1+(l0/u(1))+((l0^2)/(u(1)^2)))*(u(1)-l0))/(3*m)]';

[tarray,uarray]=ode45(f1,tspan,u0);

A=[0 -1;k0 0];
[T,U]=ode45(@(t,u)-A*u,tspan,u0);

figure(2)
plot(tarray,uarray(:,1),'LineWidth',1.2,'DisplayName',"x1 ODE45")
hold on
title("Graph for 2b)")
xlabel('t')
ylabel('l')
grid on

plot(T,U(:,1),'LineWidth',1.2,'DisplayName',"x1 ODE45")
hold on
legend('l(t)','l(t) response');

%2c)

l0=[0.05,0.1,0.2,0.5];
for ii=1:length(l0)
    f1=@(t,u) [u(2), -1*((k0)*(1+(l0(ii)/u(1))+((l0(ii)^2)/(u(1)^2)))*(u(1)-l0(ii)))/(3*m)]';

[tarray,uarray]=ode45(f1,tspan,u0);

A=[0 -1;k0 0];
[T,U]=ode45(@(t,u)-A*u,tspan,u0);

figure(ii+2)
plot(tarray,uarray(:,1),'LineWidth',1.2,'DisplayName',"x1 ODE45")
hold on
title("Graph for 2c)")
xlabel('t')
ylabel('l')
grid on

plot(T,U(:,1),'LineWidth',1.2,'DisplayName',"x1 ODE45")
hold on
legend('l(t)','l(t) response');

end