Use yyaxis
syms r
psi_list = [1, 0.1, 0.01, 0.001]; % example: 0.1, 0.01, 0.001
figure
yyaxis left
hold on
for i = 1:numel(psi_list)
psi = psi_list(i);
alpha = (12*psi)^(1/3);
l1 = @(r) (r + alpha)*(r + alpha);
l1a = @(r) r*r - r*alpha + alpha*alpha;
L1 = @(r) log(l1(r)/l1a(r));
t1 =@(r) (2*r - alpha) / (alpha*sqrt(3));
T1 =@(r) atan(t1(r));
l2 =@(r) (1 + alpha)*(1 + alpha);
l2a =@(r) (1*1) - (1*alpha) + (alpha*alpha);
L2 =@(r) log(l2(r)/l2a(r));
t2 =@(r) (2*1 - alpha) / (alpha*sqrt(3));
T2 =@(r) atan(t2(r));
W1 =@(r) 4 * sqrt(3) * (T1(r) - T2(r)) * log((r+alpha)/(1+alpha));
W2 =@(r) (L1(r) - L2(r))*log((r*r - r*alpha + alpha*alpha)/(1 - alpha + alpha*alpha));
W3 =@(r) 4*(T1(r) -T2(r))*(T1(r) -T2(r));
W4 =@(r) (L1(r) - L2(r)) + 2 * sqrt(3) * (T1(r) - T2(r));
W5 =@(r) 3 ./ ((r - 1) * (r - 1));
W =@(r) ((W1(r) + W2(r) -W3(r)) * W5(r))/W4(r);
F =@(r) ((W(r) * (r - 1) * (r - 1) / 2) + log(r) ) / (r - 1);
ylim([0.5 1])
xlim([1 6])
set(gca, 'ylim', [0.5 1]);
set(gca, 'ytick', 0.5:0.1:1);
set(gca, 'xlim', [1 6]);
set(gca, 'xtick', 1:1:6);
fplot(F(r), [1 6])
end
yyaxis right
hold on
for i = 1:numel(psi_list)
psi = psi_list(i);
alpha = (12*psi)^(1/3);
l1 = @(r) (r + alpha)*(r + alpha);
l1a = @(r) r*r - r*alpha + alpha*alpha;
L1 = @(r) log(l1(r)/l1a(r));
t1 =@(r) (2*r - alpha) / (alpha*sqrt(3));
T1 =@(r) atan(t1(r));
l2 =@(r) (1 + alpha)*(1 + alpha);
l2a =@(r) (1*1) - (1*alpha) + (alpha*alpha);
L2 =@(r) log(l2(r)/l2a(r));
t2 =@(r) (2*1 - alpha) / (alpha*sqrt(3));
T2 =@(r) atan(t2(r));
W1 =@(r) 4 * sqrt(3) * (T1(r) - T2(r)) * log((r+alpha)/(1+alpha));
W2 =@(r) (L1(r) - L2(r))*log((r*r - r*alpha + alpha*alpha)/(1 - alpha + alpha*alpha));
W3 =@(r) 4*(T1(r) -T2(r))*(T1(r) -T2(r));
W4 =@(r) (L1(r) - L2(r)) + 2 * sqrt(3) * (T1(r) - T2(r));
W5 =@(r) 3 ./ ((r - 1) * (r - 1));
W =@(r) ((W1(r) + W2(r) -W3(r)) * W5(r))/W4(r);
fplot(W(r), [1 6])
end
hold off
ylim([0 0.2])
xlim([1 6])
set(gca, 'ylim', [0 0.2]);
set(gca, 'ytick', 0:0.05:0.2);
set(gca, 'xlim', [1 6]);
set(gca, 'xtick', 1:1:6);
