Why the following code is not generating plot for W?

Reference figure: plot should be like this for W refer square ones.

Ho = 1;
alpha = 0.1;

a = 1.0;
sigma = 0.15;

lbar_list = [0.0, 0.2];
eps = 0:0.2;
hold on
for i = 1:numel(lbar_list)
    lbar = lbar_list(i);
    H = @(x) Ho + a*(1 - x);            
    G1 = @(x,eps) H(x).^3 + 3 .* H(x).^2 .* alpha + 3 .* H(x) .* alpha^2 + 3 .* H(x) .* sigma^2 + eps + 3*sigma^2*alpha + alpha^3 - 12.*lbar^2 .* (H(x) + alpha);
    G2 = @(x) 24 .* lbar.^3 .* tanh(H(x)./(2.*lbar));
    G3 = @(x,eps) (12.*lbar^2*alpha - eps - alpha^3 - 3*sigma^2*alpha) .* (1 - (tanh(H(x)./(2.*lbar))).^2);
    G = @(x,eps) G1(x,eps) + G2(x) + G3(x,eps);
    Hm1 = @(x,eps) H(x).* (1 ./ G(x,eps));
    Hm2 = @(x,eps) (1 ./ G(x,eps));
    IntHm1 =@(eps) integral(@(x) Hm1(x,eps),0,1);
    IntHm2 =@(eps) integral(@(x) Hm2(x,eps),0,1);
    Hm = @(eps) IntHm1(eps) / IntHm2(eps);
    P1 = @(x,eps) 6 .* (1 ./ G(x,eps)) .* (H(x) - Hm(eps));
    P2 = @(x,eps) integral(@(x) P1(x,eps),0,x);
   
    W = @(eps) integral(@(x) P2(x,eps),0,1, 'ArrayValued', true);

    plot(eps,4.5*W(eps))

end

legend(num2str(lbar_list'))
ylim([0.4 0.8])
xlim([0 0.2])
set(gca, 'ytick', 0.4:0.1:0.8);
set(gca, 'xtick', 0:0.02:0.2);
xlabel('eps')
ylabel('W')

Answer 1

[KSSV]

I remember you asking the similiar question earlier. What have you learned from the previous question? You are keep on asking same questions. Spend some time on the documentation and rethink on your old codes.

clc; clear all ;

Ho = 1;
alpha = 0.1;

a = 1.0;
sigma = 0.15;

lbar_list = [0.0, 0.2];
eps = linspace(0,0.2);
hold on
for i = 1:numel(lbar_list)
    lbar = lbar_list(i);
    H = @(x) Ho + a*(1 - x);            
    G1 = @(x,eps) H(x).^3 + 3 .* H(x).^2 .* alpha + 3 .* H(x) .* alpha^2 + 3 .* H(x) .* sigma^2 + eps + 3*sigma^2*alpha + alpha^3 - 12.*lbar^2 .* (H(x) + alpha);
    G2 = @(x) 24 .* lbar.^3 .* tanh(H(x)./(2.*lbar));
    G3 = @(x,eps) (12.*lbar^2*alpha - eps - alpha^3 - 3*sigma^2*alpha) .* (1 - (tanh(H(x)./(2.*lbar))).^2);
    G = @(x,eps) G1(x,eps) + G2(x) + G3(x,eps);
    Hm1 = @(x,eps) H(x).* (1 ./ G(x,eps));
    Hm2 = @(x,eps) (1 ./ G(x,eps));
    IntHm1 =@(eps) integral(@(x) Hm1(x,eps),0,1);
    IntHm2 =@(eps) integral(@(x) Hm2(x,eps),0,1);
    Hm = @(eps) IntHm1(eps) / IntHm2(eps);
    P1 = @(x,eps) 6 .* (1 ./ G(x,eps)) .* (H(x) - Hm(eps));
    P2 = @(x,eps) integral(@(x) P1(x,eps),0,x);
   
    W = @(eps) integral(@(x) P2(x,eps),0,1, 'ArrayValued', true);
    WW = zeros(1,length(eps)) ; 
    for j = 1:length(eps)
       WW(j) = 4.5*W(eps(j)) ; 
    end
    plot(eps,WW)

end

legend(num2str(lbar_list'))
ylim([0.4 0.8])
xlim([0 0.2])
set(gca, 'ytick', 0.4:0.1:0.8);
set(gca, 'xtick', 0:0.02:0.2);
xlabel('eps')
ylabel('W')

[Torsten]

@AVINASH SAHU comment moved here:

Okay thanks I will spend more time reading the documentation