L=80; %Sites

%s=1/8;
%lambda=s*L*L;

syms H K a u;

a=(2*sqrt(lambda*(2*H+lambda)))/(H+lambda+sqrt(lambda*(2*H+lambda)));
u=1-a;

%K=2*log(2)-log(u)/2+u*(-1*log(u)-2+4*log(2))/8;
K=ellipke(a);
eqn = 1/2 == ((sqrt(2*H+lambda)-sqrt(lambda))*K)/(sqrt(2)*H);
cond = (lambda+sqrt(lambda*(2*H+lambda)))/(2*lambda) < 1;
assume(cond);
H=solve(eqn,H); %Calculate Hamiltonian

b=(2*sqrt(lambda*(2*H+lambda)))/(H+lambda+sqrt(lambda*(2*H+lambda)));
v=1-b;
%%%E=1+v*(-1*log(v)/4-1/4+log(2))+((v^2)/64)*(-6*log(v)-13+24*log(2));
[K,E]=ellipke(b);

phi=-1*sqrt(2)*(sqrt(2*H+lambda)+sqrt(lambda))*E;
psi=double(phi/L);
philist(i+1)=phi;
lamlist(i+1)=20+(i/100);