%ONESIM_BIGGER.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.
% 
% Uses a larger system size, N, than onesim.m.
%
% DJH Sep 2019

r = 0.5;
theta = 0.5;

N = 100;

% 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 = 200-1;
Usol = zeros(N,steps+1);
Usol(:,1) = Uzero;
Udiff = zeros(steps,1); % record change in solution over each time step
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;
   Udiff(k) = norm(Ukp1-Uk)/norm(Uk);
   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_biggersima.pdf 
print -dpng pic_biggersima.png

figure(2)
clf
semilogy(Udiff,'LineWidth',4)
xlim([1 steps])
xlabel('Time')
ylabel('Relative change in U') 
set(gca,'FontSize',14,'FontWeight','Bold')
print -dpdf pic_biggersimb.pdf 
print -dpdf pic_biggersimb.png 

figure(3)
clf
plot(Usol(:,end),'k-v','LineWidth',4,'MarkerSize',10)
xlim([1 N])
xlabel('U component')
ylabel('opinion')
legend('Final','Location','Northwest') 
set(gca,'FontSize',14,'FontWeight','Bold')
print -dpdf pic_biggersimc.pdf 
print -dpdf pic_biggersimc.png 