%ONESIM.M
% Accompanies BIFSIM.M
% DOES A SINGLE RUN WHERE U1 and U2 ARE STABLE.
%
% Shows how the existence of two different, stable solutions can 
% affect the intermediate time behaviour.
% Eventually, one of the two solutions dominates. 
% In this case the lower value, u1, wins out in the end.
% But we can see u2 trying to attract the solution before it loses out.
%
% DJH Dec 2016
%
% This version Oct 2017
% Now Sep 2019

r = 0.5;
theta = 0.5;

N = 32;

% Set up A for nearest neighbour ring network
A = zeros(N,N);
A = diag(ones(N-1,1),1) + diag(ones(N-1,1),-1);
A(1,N) = 1;
A(N,1) = 1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

epsi = 2.5;
u1 = 0.5 - 0.5*sqrt(1 - 4*r/epsi)
u2 = 0.5 + 0.5*sqrt(1 - 4*r/epsi)

%%%%% Initial values %%%%
x = linspace(-10,10,N)'; 
sig1 = exp(x)./(1 + exp(x));
Uzero = 0.3 + 0.225*sig1;
Uk = Uzero;
steps = 500-1;
steps = 300-1;
Usol = zeros(N,steps+1);
Usol(:,1) = Uzero;
for k = 1:steps   
   Ukp1 = Uk + r*(theta - Uk) + (epsi)*(Uk.*(1-Uk)).*(0.5*A*Uk - theta);      
             %%%%%% Reflection---could make this better? %%%%
   [a1,b1] = find(Ukp1<0);
   [a2,b2] = find(Ukp1>1);
   Ukp1(a1) = 0;        % just perturb back to 0 
   Ukp1(a2) = 1;        % just perturb back to 1
   Usol(:,k+1) = Ukp1;
   Uk = Ukp1;
end

figure(1)
clf
waterfall(Usol')
ylabel('time')
xlabel('U component')
zlabel('opinion')
xlim([1 N])
set(gca,'FontSize',14,'FontWeight','Bold')
print -dpdf pic_onesima.pdf 
print -dpng pic_onesima.png

figure(2)
clf
plot(Usol(:,1),'r-o','LineWidth',4,'MarkerSize',10)
hold on
plot(Usol(:,(steps+1)/2),'b-^','LineWidth',4,'MarkerSize',10)
plot(Usol(:,end),'k-v','LineWidth',4,'MarkerSize',10)
xlim([1 N])
xlabel('U component')
ylabel('opinion')
legend('Initial','Intermediate','Final','Location','Northwest') 
set(gca,'FontSize',14,'FontWeight','Bold')
print -dpdf pic_onesimb.pdf 
print -dpdf pic_onesimb.png 
