How to plot these two plots in a single figure and using both y axis?

syms r
psi_list = [1, 0.1, 0.01, 0.001]; % example: 0.1, 0.01, 0.001
figure
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

figure
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);

Answer 1

[Chunru]

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);

[AVINASH SAHU]

thank you so much