Write a matlab code which calculates and plots the change of the potential inside and outside of the spherical Earth. Mass of the Earth and its radius are known. Consider Earth interior is homogeneous with the same density

Okay so ive tried working on the code and this is what i got thus far. when i run my code, it plots a figure with no curve. I am supposed to get a graph that looks like this. Please help.

%This code will calculate and plots the change in gravity pontential inside
%and outside of the spherical Earth

%Let U represent gravitational potential
%G = gravitational constant (m^3kg^-1s^-2)
G = 6.67*10E-11;
%M = mass of Earth (kg)
M = 5.97*10E24;
%r = distance from center of Earth to point mass (m)
%Let's assume a point mass is located at a radius of r = 1000000 m and
%plot the results
r = 10000000;
%R = radius to the surface of the Earth (m)
R = 6371000;
%p = density (kg/m^3)

%Using Poisson's equation we can derive a formula that can give us the
%gravitational potential inside the spherical Earth. Let U1 represent the
%potential inside the sphere.
for x = 0:1000000:r 
    while r < R
        U1 = -(G*M)*(3*R^2-r^2)/(2*R^3);
        break
    end
end

%We can derive Laplace's equation to find an equation that will give us the
%gravitational potential at a point on the surface of the Earth. Let U2 
%represent the potential on the Earth.
for x = 0:1000000:r 
    while r > R
        U2 = -(G*M)/r; 
        break
    end;
end;

%The potential at the center of sphere is obtained by putting r = 0: Let Uc
%represent the potential at the centre of the Earth.

Uc = -(3*G*M)/2*r;
v = r/R;
plot(v,U1,v,U2,v,Uc);
xlabel('r/R');
ylabel('U(m^2s^2)');
title('Plot of Gravity Potential as a function of radius');
set(gca,'XAxisLocation','top');

[Walter Roberson]

When you edit away your question, the volunteers usually stop answering you. This is a shared resource: we answer for everyone to learn. If you want to ask your questions privately you can hire someone to answer them.

[Rena Berman]

(Answers Dev) Restored question.

Answer 1

[Walter Roberson]

Look at your code

for x = 0:1000000:r 
    while r < R
        U1 = -(G*M)*(3*R^2-r^2)/(2*R^3);
        break
    end
end

Your loop is a "for" loop over x, but you do not use x inside the body.

Your initial r > R so the while loop exits immediately every time.

If your initial r < R then you are writing a value of U1 for each iteration but the value does not depend upon the previous values, and you do not save the U1 values as you go. The effect is as-if you wasted time for a while and then only ran the last iteration of the loop.

You have similar problems for your second "for".

When you get to the end to calculate v in terms of r and R, your r and R only have a single value, so your v is going to be scalar.

[vimal kumar chawda]

Can you share the right code as I am facing the same problem? For the same task please.