function cella_benard(x0,Ac,B,A)
% programma simulazione cella Benard
% x0=1;  velocità iniziale del fluido
% B=1;
% Ac=1;
% A=0.5; poichè (A<Ac) -> un pto di eq., x=0. 
% qui DT<valore di soglia cioè diventa immonbile
% A=2; poichè (A>Ac) -> 3 pti di eq., x=0,  x=-1, x=+1
% qui DT>valore di soglia cioè liquido rimane in  moto
% punti di equilibrio se A = 2:
% x1 = 0; x2 = +/-sqrt((A-Ac)/B)) = +/-sqrt((2-1)/1))=+/-1

if nargin < 4, 
    fprintf(' inserire x0, Ac, B, A\n')
    return
end

if A==Ac
    time=5*B;
else
    time=2*abs(1/(A-Ac));
end
DT=time/10000;
N=floor(time/DT);
x=zeros(N,1);
x(1)=x0;
for i = 1:N
    x1=((A-Ac)*x(i)-B*x(i)^3)*DT;
    xn=x(i)+x1/2;
    x2=((A-Ac)*xn-B*xn^3)*DT;
    xn=x(i)+x2/2;
    x3=((A-Ac)*xn-B*xn^3)*DT;
    xn=x(i)+x3;
    x4=((A-Ac)*xn-B*xn^3)*DT;
    x(i+1,1)=x(i)+x1/6+x2/3+x3/3+x4/6;
end

if A-Ac > 0,
    peq=sqrt((A-Ac)/B);
    xmax=peq*1.5;
else
    peq = 0;
    xmax=max(abs(x))*1.5;
end

if abs(xmax)<abs(x0)
    xmax=x0;
end

t=(0:DT:time);
figure(1);
plot(t, x);
ylabel('Velocità fluido')
xlabel('tempo');
s_title=sprintf('Cella Benard RK4');
title(s_title);
%xlim([t(1) t(end)]);
%ylim([peq max(x)*1.2]);

fprintf('2° punto di equilibrio: +/-%.3e\n', peq)
DX=xmax/1000;
x=(-xmax:DX:xmax)';
U=-(A-Ac)/2*x.^2+B/4*x.^4;
U0=-(A-Ac)/2*x0^2+B/4*x0^4;
U1=-(A-Ac)/2*peq^2+B/4*peq^4;
figure(2)
  plot(x,U, x0, U0, 'ro', 0,0, 'rx', peq, U1, 'rx', -peq, U1, 'rx')
  ylabel('Potenziale')
  xlabel('Velocità');
  s_title=sprintf('Benard potenziale efficace');
  title(s_title);

return