Please find the mistake in this code.

1 view (last 30 days)
Reference image: plot should be like this but don't know where I am making mistake.
syms x lambda
warning off
alpha = -0.1;
sigma = 0.1;
eps = -0.1;
e = 0.2;
M = 10;
a = 2;
figure
psi_list = [10, 20, 30, 40. 50];
for i = 1:numel(psi_list)
psi = psi_list(i);
hbar = @(x) a - a.*x + x;
A1 = eps + alpha^3 + (3 * sigma^2 * alpha);
B1 =@(x,lambda) (-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 = @(lambda) 2 * (M^2) * (1 + lambda);
C1=@(x) a1(x) + (b1(x) .* c1);
D1 =@(x,lambda) d1(lambda) .* ((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,lambda) B1(x,lambda) + (3 * lambda .* C1(x) .* hbar(x)) + (3 * lambda .* C1(x) .* alpha) + (D1(x,lambda) .* C1(x));
f2 =@(x,lambda) 12 * (M^2) * (1 + lambda) .* C1(x);
f3 = psi * (e^3);
f4 = @(lambda) (1 + lambda) *180 * ((1 - e)^2);
f5 =@(lambda) 1/(2 + lambda);
F = @(x,lambda) ((f5(lambda) .* f1(x,lambda))./f2(x,lambda)) + (f3/f4(lambda));
q1 =@(x,lambda) hbar(x) ./ (2 .* F(x,lambda));
Q1 = int(q1,x,0,1); % integration with respect to "x" limit 0 to 1
q2 = @(x,lambda) 1./(F(x,lambda));
Q2 = int(q2,x,0,1); % integration with respect to "x" limit 0 to 1
Q = @(lambda) Q1(lambda)/Q2(lambda); % Q is a function of "lambda"
p1 = @(x,lambda) (1./F(x,lambda)) .* ((0.5 .* hbar(x)) - Q(lambda));
P = @(x,lambda) int(p1,x); % integration with respect to "x" limit 0 to x
% W is obtained by integrating "P" with respect to x limit to x is 0 to
% 1
W =@(lambda) int(P,x,0,1);
fplot(@(lambda) W(lambda))
end
legend(num2str(psi_list'))
ylim([0 1])
xlim([0 4])
set(gca, 'ytick', 0:0.1:1);
set(gca, 'xtick', 0:0.2:4);
xlabel('lambda')
ylabel('W(lambda)')
  6 Comments

Sign in to comment.

Accepted Answer

Torsten
Torsten on 8 Jul 2022
Edited: Torsten on 8 Jul 2022
alpha = -0.1;
sigma = 0.1;
eps = -0.1;
e = 0.2;
M = 10;
a = 2;
%figure
psi_list = [10, 20, 30, 40, 50];
lambda = 0:0.2:4;
for i = 1:numel(psi_list)
psi = psi_list(i);
hbar = @(x) a - a.*x + x;
A1 = eps + alpha^3 + (3 * sigma^2 * alpha);
B1 =@(x,lambda) (-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 = @(lambda) 2 * (M^2) * (1 + lambda);
C1=@(x) a1(x) + (b1(x) .* c1);
D1 =@(x,lambda) d1(lambda) .* ((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,lambda) B1(x,lambda) + (3 * lambda .* C1(x) .* hbar(x)) + (3 * lambda .* C1(x) .* alpha) + (D1(x,lambda) .* C1(x));
f2 =@(x,lambda) 12 * (M^2) * (1 + lambda) .* C1(x);
f3 = psi * (e^3);
f4 = @(lambda) (1 + lambda) *180 * ((1 - e)^2);
f5 =@(lambda) 1./(2 + lambda);
F = @(x,lambda) ((f5(lambda) .* f1(x,lambda))./f2(x,lambda)) + (f3./f4(lambda));
q1 =@(x,lambda) hbar(x) ./ (2 .* F(x,lambda));
Q1 = @(lambda)integral(@(x)q1(x,lambda),0,1,'ArrayValued',true); % integration with respect to "x" limit 0 to 1
q2 = @(x,lambda) 1./(F(x,lambda));
Q2 = @(lambda)integral(@(x)q2(x,lambda),0,1,'ArrayValued',true); % integration with respect to "x" limit 0 to 1
Q = @(lambda) Q1(lambda)./Q2(lambda); % Q is a function of "lambda"
p1 = @(x,lambda) (1./F(x,lambda)) .* ((0.5 .* hbar(x)) - Q(lambda));
P = @(x,lambda) integral(@(y)p1(y,lambda),0,x,'ArrayValued',true); % integration with respect to "x" limit 0 to x
% W is obtained by integrating "P" with respect to x limit to x is 0 to
% 1
W =@(lambda) integral(@(x)P(x,lambda),0,1,'ArrayValued',true);
%fplot(@(lambda) W(lambda))
plot(lambda,W(lambda))
hold on
end
legend(num2str(psi_list'))
ylim([0 1])
xlim([0 4])
set(gca, 'ytick', 0:0.1:1);
set(gca, 'xtick', 0:0.2:4);
xlabel('lambda')
ylabel('W(lambda)')
  1 Comment
AVINASH SAHU
AVINASH SAHU on 8 Jul 2022
Thank you so much. It helped me a lot

Sign in to comment.

More Answers (0)

Categories

Find more on Creating and Concatenating Matrices in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!