r= 3.2
 s= 0.03
 f= 0.8
%% Find stationary states
syms xi
eqn = -(r^2 * (xi-1)*((xi/2)-1)+1+(s/(1-xi)))*xi + f ;
%eqn2= -(r^2 * (xi-1)*((xi/2)-1)*(1-xi)+(1-xi)+(s))*xi+(f*(1-xi));

f1 = matlabFunction(expand(eqn))

root1=fzero(f1,[-4, 0.999999999])
root2=fzero(f1,[1.0000001, 4])


%% Solve ODE
omega = 0.01;
T = pi/omega;

t_start = 0;
t_end = T;
time_span = t_start:0.01:t_end;

xi0 = [0 0];

[t,xi] = ode45(@odefun, time_span,xi0);
plot(t,xi);
hold on
plot (t, 3*sin(omega*t));
ylim([0,3.5])
legend('Displacement','Force','location','NorthEast')
title('Displacement vs Time')
xlabel('Time, tau [sec]')
ylabel('Displacement, \xi')
yyaxis right
ylabel ('Force F/F0')
ylim([0,3.5])
function xidot = odefun(t,xi)
xidot = -((3.2^2).*(xi-1).*(0.5*xi-1).*xi+xi+(xi.*0.03)./(1-xi))+3*sin(0.01*t);
end