Can someone help me check if I did this correctly?

Here's the problem:

And the script:

syms r a0;

% Part (a) - Probability Density Function p(r) and its integral from 0 to infinity
p_r = 4 / a0^3 * r^2 * exp(-2 * r / a0);
integral_p = int(p_r, r, 0, Inf);
disp("Part (a) - Integral of p(r) from 0 to infinity:");

Part (a) - Integral of p(r) from 0 to infinity:

disp(integral_p);

% Part (b) - Limit as r approaches infinity
limit_r_infinity = limit(p_r * r^2 / exp(2 * r / a0), r, Inf);
disp("Part (b) - Limit as r approaches infinity:");

Part (b) - Limit as r approaches infinity:

disp(limit_r_infinity);

% Part (c) - Visualization of p(r) and its maximum
a0 = 5.59e-11;
p_r = @(r) 4 / a0^3 * r.^2 .* exp(-2 * r / a0);
r_vals = linspace(0, 2*a0, 1000);
p_values = p_r(r_vals);
[max_value, max_index] = max(p_values);
max_r = r_vals(max_index);

plot(r_vals, p_values);
hold on;
plot(max_r, max_value, 'ro');
text(max_r, max_value, sprintf('(%.2e, %.2e)', max_r, max_value), 'VerticalAlignment', 'bottom');
title('Probability Density Function p(r)');
xlabel('r');
ylabel('p(r)');
legend('p(r)', 'Maximum');
hold off;

% Part (d) - Calculation of P(4a0)
P_4a0 = integral(p_r, 0, 4 * a0);
disp("Part (d) - Integral P(4a0):");

Part (d) - Integral P(4a0):

disp(P_4a0);

    0.9862

% Part (e) - Expected value mu
syms r a0;

% Define the integrand for the expected value mu
integrand_mu = r^3 * exp(-2 * r / a0);

% Calculate the expected value mu
mu = 4 / a0^3 * int(integrand_mu, r, 0, Inf);
disp("Part (e) - Expected value mu:");

Part (e) - Expected value mu:

disp(mu);

Answer 1

[Torsten]

Replace

syms r a0 

by

syms r a0 positive

to get the limits resolved.