hold on
colorstring = 'kbgrc';
t = linspace(0,0.617,1000000);

% In the code, B corresponds to 1/gamma_c and A corresponds to gamma_f
% simulate progressively for different endtimes in increasing order without
% closing the figure window
% endtimes - 0.070, 0.283,0.779, 2.444 for B=100 
% endtimes - 0.069, 0.264,0.5859, 0.617 for B=1
% endtimes  - 0.0061, 0.0078, 0.0102, 0.0138 for B=0.01
B=1;
%rho_0 =  B+1-sqrt(B^2-2*B*t)
%plot(t,rho_0,'Color','y')

for i =1:4 %(if B=0.01, this should be progressively changed from 1:4 to 2:4,3:4 and 4:4, as the endtime is increased)
    
A = 0.2+(i-1)*0.2
%nu = 1+i/2; %change this
%line([2*B/((2-A)^2) 2*B/((2-A)^2)], [1 1.06], 'Color',colorstring(i), 'LineStyle','--'); % this is the right option for B=0.01
line([2*B/((2-A)^2) 2*B/((2-A)^2)], [1 5], 'Color',colorstring(i), 'LineStyle','--');
q = 1 - A/2;
F1 = q*(1+A*B*q*t).^2+ B*(1-A*B)^2*q*t + (B-A*B*B)*(1+A*B*q*t);
F2 = 1+A*B*q*t;
F3 = B*(1-A*B)^2*(2*q*t+2*A*B*q*q*t.^2);
F4 = B^2*(1-A*B)^2*(1+q*t).^2;
%(2*F1.*F2-F3).^2 + 4*F2.^2.*(F4-F1.^2)
F5 = (F3-2*F1.*F2 + sign(1-A*B)* realsqrt(F3.^2 + 4*F4.*F2.^2 -4*F1.*F2.*F3))./(2*F2.^2);
rho = 1/A*(1-realsqrt(1+2*A*F5));

%
plot(t,rho, 'Color', colorstring(i),'Linewidth',2)
line([0 2.44], [1/A 1/A], 'Color',colorstring(i), 'LineStyle','--'); %comment this out for B=0.01
xlim([0 0.8])  % this should range up to 0.015 for B=0.01, 0.8 for B=1, 3 for B=100
end