I'm very upset this code doen't run a result,it just running all the time.I think the reason is the analytic derivation process which need a analytic solution of matrix W. I hope someone would give me some help or point the problems,I will be very grateful.
tic;
clc;
clear;
syms B
W11=fa(B,1,1,1,1);W12=fa(B,1,1,1,0);W13=fa(B,1,1,0,1);W14=fa(B,1,1,0,0);
W21=fa(B,1,0,1,1);W22=fa(B,1,0,1,0);W23=fa(B,1,0,0,1);W24=fa(B,1,0,0,0);
W31=fa(B,0,1,1,1);W32=fa(B,0,1,1,0);W33=fa(B,0,1,0,1);W34=fa(B,0,1,0,0);
W41=fa(B,0,0,1,1);W42=fa(B,0,0,1,0);W43=fa(B,0,0,0,1);W44=fa(B,0,0,0,0);
W=simplify([W11,W12,W13,W14;W21,W22,W23,W24;W31,W32,W33,W34,W41,W42,W43,W44]);
E1=eig(W);
E2=sort(E1);
E3=E2(4);
T=1;
be=1/T;
f1=-1/be*log(E3);
ma=-diff(f1);
m1=subs(ma,{B},2)
toc;
function s1=fa(B,m1i,m1j,m2i,m2j)
JH=2;J=1;T=1;de=0.2;
sx=1/2*[0 1;1 0];
sy=1/2*[0 -1i;1i 0];
sz=1/2*[1 0;0 -1];
u=eye(2);
be=1/T;
H11=kron(sx,kron(sx,kron(u,u)))+kron(sy,kron(sy,kron(u,u)))...
+de*kron(sz,kron(sz,kron(u,u)));
H12=kron(u,kron(sx,kron(sx,u)))+kron(u,kron(sy,kron(sy,u)))...
+de*kron(u,kron(sz,kron(sz,u)));
H2=kron(sz,kron(u,kron(u,u)))+kron(u,kron(sz,kron(u,u)));
H3=kron(u,kron(u,kron(u,u)));
H=-JH*(H11+H12)+J*H2*(m1i+m1j)+1/2*J*H3*(m1i+m2i+m1j+m2j)...
-B*1/2*H3*(m1i+m2i+m1j+m2j);
va=eig(H);
s1=sum(exp(-be*va));
end
As suggested by @Dev-iL in a comment, the problem with the code is that it tries to compute very complicated expressions symbolically. It is much simpler to evaluate the expression numerically instead. The numerical approximation to the derivative will be the only place where you don't get exact results, but as we'll see later, the error is very, very small.
This is how I rewrote your code. First I took all your code up to the diff line, and made it into a function (simplified slightly for clarity). This function does not use symbolic math at all, it simply computes the value of f1 at the given value of B. Note that it uses the fa function in the question unchanged.
function out = f1(B)
W = [fa(B,1,1,1,1), fa(B,1,1,1,0), fa(B,1,1,0,1), fa(B,1,1,0,0);...
fa(B,1,0,1,1), fa(B,1,0,1,0), fa(B,1,0,0,1), fa(B,1,0,0,0);...
fa(B,0,1,1,1), fa(B,0,1,1,0), fa(B,0,1,0,1), fa(B,0,1,0,0);...
fa(B,0,0,1,1), fa(B,0,0,1,0), fa(B,0,0,0,1), fa(B,0,0,0,0)];
E = sort(eig(W));
E = E(4);
T = 1;
be = 1/T;
out = -1/be*log(E);
end
Next, we estimate the derivative numerically by computing the value of f1 at two locations really close to the point B (but not too close, as we'll get into numerical rounding errors):
B = 2;
delta = 1e-6;
(f1(B-delta) - f1(B+delta))/(2*delta)
We change the value of delta to verify we have a good approximation. With 1e-6 and 1e-4 and 1e-3 I get the exact same value: 1.6303. This indicates that the function is very smooth around B, and our estimate of the derivative is correct.