Why it is not making plot for different values of "M"?

syms x
alpha = -0.1;
sigma = 0.1;
eps = -0.1;
e = 0.2;
a = 2;
lambda = 2;

psi = 10;

figure
M_list = [4, 6, 8, 10];
for i = 1:numel(M_list)
    M = M_list(i);
    hbar = @(x) a - a.*x  + x;
    A1 = eps + alpha^3 + (3 * sigma^2 * alpha);
    B1 = @(x) (-3 * lambda * M) * ((hbar(x).^2) + (2 .* hbar(x) .* alpha) + (sigma^2) + (alpha^2));

    a1 = @(x) tanh(M .* hbar(x));
    b1 = @(x) 1 - ((tanh(M .* hbar(x))).^2);
    c1 = (M * alpha) - ((M^3 * A1)/3);
    d1 = 2 * (M^2) * (1 + lambda);  
    
    C1 = @(x) a1(x) + (b1(x) .* c1);
    D1 = @(x) d1 .* ((hbar(x).^3) + (3 .* (hbar(x).^2) .* alpha) + (3 .* hbar(x) .* (alpha)^2) + (3 .* hbar(x) .* (sigma)^2) + eps + (3 * alpha * (sigma^2)) + (alpha^3));

    f1 = @(x) B1(x) + (3 * lambda .* C1(x) .* hbar(x)) + (3 * lambda .* C1(x) .* alpha) + (D1(x) .* C1(x));
    f2 = @(x) 12 * (M^2) * (1 + lambda) .* C1(x);
    f3 = psi * (e^3);
    f4 = (1 + lambda) *180 * ((1 - e)^2);  % 180 is not given in paper 
    f5 = 1/(2 + lambda);

    F = @(x) ((f5 .* f1(x))./f2(x)) + (f3/f4);

    q1 = @(x) hbar(x) ./ (2 .* F(x));
    Q1 = integral(q1,0,1);

    q2 = @(x) 1./(F(x));
    Q2 = integral(q2,0,1);

    Q = Q1/Q2;

    p1 = @(x) (1./F(x)) .* ((0.5 .* hbar(x)) - Q);
    P = @(x) integral(p1,0,x);

    fplot(P, [0 1])
    ylim([0 1])
    set(gca, 'ytick', 0:0.1:1);
    set(gca, 'xtick', 0:0.2:1);
    xlabel('x')
    ylabel('P(x)')
    
end

Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is 1.6e+02. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.

Warning: Function behaves unexpectedly on array inputs. To improve performance, properly vectorize your function to return an output with the same size and shape as the input arguments.

Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is 8.0e+02. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.

Warning: Function behaves unexpectedly on array inputs. To improve performance, properly vectorize your function to return an output with the same size and shape as the input arguments.

Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is 1.7e+02. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.

Warning: Function behaves unexpectedly on array inputs. To improve performance, properly vectorize your function to return an output with the same size and shape as the input arguments.

Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is 2.1e+02. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.

Warning: Function behaves unexpectedly on array inputs. To improve performance, properly vectorize your function to return an output with the same size and shape as the input arguments.

Answer 1

[KSSV]

You have to use hold on.

syms x
warning off
alpha = -0.1;
sigma = 0.1;
eps = -0.1;
e = 0.2;
a = 2;
lambda = 2;

psi = 10;


M_list = [4, 6, 8, 10];
figure
hold on
for i = 1:numel(M_list)
    M = M_list(i);
    hbar = @(x) a - a.*x  + x;
    A1 = eps + alpha^3 + (3 * sigma^2 * alpha);
    B1 = @(x) (-3 * lambda * M) * ((hbar(x).^2) + (2 .* hbar(x) .* alpha) + (sigma^2) + (alpha^2));
    
    a1 = @(x) tanh(M .* hbar(x));
    b1 = @(x) 1 - ((tanh(M .* hbar(x))).^2);
    c1 = (M * alpha) - ((M^3 * A1)/3);
    d1 = 2 * (M^2) * (1 + lambda);
    
    C1 = @(x) a1(x) + (b1(x) .* c1);
    D1 = @(x) d1 .* ((hbar(x).^3) + (3 .* (hbar(x).^2) .* alpha) + (3 .* hbar(x) .* (alpha)^2) + (3 .* hbar(x) .* (sigma)^2) + eps + (3 * alpha * (sigma^2)) + (alpha^3));
    
    f1 = @(x) B1(x) + (3 * lambda .* C1(x) .* hbar(x)) + (3 * lambda .* C1(x) .* alpha) + (D1(x) .* C1(x));
    f2 = @(x) 12 * (M^2) * (1 + lambda) .* C1(x);
    f3 = psi * (e^3);
    f4 = (1 + lambda) *180 * ((1 - e)^2);  % 180 is not given in paper
    f5 = 1/(2 + lambda);
    
    F = @(x) ((f5 .* f1(x))./f2(x)) + (f3/f4);
    
    q1 = @(x) hbar(x) ./ (2 .* F(x));
    Q1 = integral(q1,0,1);
    
    q2 = @(x) 1./(F(x));
    Q2 = integral(q2,0,1);
    
    Q = Q1/Q2;
    
    p1 = @(x) (1./F(x)) .* ((0.5 .* hbar(x)) - Q);
    P = @(x) integral(p1,0,x);
    
    fplot(P, [0 1])
    
    
end
legend(num2str(M_list'))
ylim([0 1])
set(gca, 'ytick', 0:0.1:1);
set(gca, 'xtick', 0:0.2:1);
xlabel('x')
ylabel('P(x)')

[AVINASH SAHU]

Thank you so much