How to graph my calculations

I am having so much trouble figuring out what is wrong with my code. The graph should look like this, but it is not working out. I believe I calculated everything correctly. density and ptemp are previous functions created. Please help!

Write a program, stability.m, to compute and plot N (buoyancy frequency) as a function of depth or pressure using the above equation (N = [(-g/RHO)(dRHO_theta/dz)]^1/2. The derivative must be computed numerically; this may take you some time to figure out if you have not done it before, so be patient. Computing derivatives numerically is a key skill in many sciences. Also compute N by substituting r and rt for rq. Do not plot any imaginary values of N you may find (replace with NaN if needed). Report the buoyancy period at depths of approximately 60, 500, and 1000 m, which will give you an idea of the period of internal waves propagating around the ocean at these depths

% using hydrostatic equation to compute z  (dp/dz = -RHO*g)
 
clear
[p t s x1 x2 x3 x4 x5] = textread('hotctddata.txt','%f %f %f %f %f %f %f %f',...
     'headerlines',6);
n = size(t,1);
 
for i = 1:n
    RHO(i) = density(t(i),s(i),p(i));
    RHOp(i) = density(t(1),s(1),p(i));
    RHOt(i) = density(t(i),s(1),p(1));
    RHOs(i) = density(t(1),s(i),p(1));
end
 
z = zeros(n,1);
for i = 2:n
    g = -9.8;
    dp = (p(i)-p(i-1))*10000;
    RHOavg(i) = (RHO(i)+RHO(i-1))*0.5;
    dz = dp./(RHOavg(i)*g);
    z(i) = dz+z(i-1);
end
 
% solving for N = sqrt((-g/RHO)*(dptheta/dz))
 
for i=1:n
  theta(i) = ptemp(t(i),s(i),p(i));
end
 
for i = 1:n
  RHO(i) = density(t(i),s(i),p(i));
  RHOtheta(i) = density(theta(i),s(i),p(1));
  RHOt(i) = density(t(i),s(i),p(1));
end
 
for i = 2:n
    g = -9.8;
    RHOavg(i) = (RHO(i)+RHO(i-1))*0.5;
    dptheta = (RHOtheta(i) - RHOtheta(i-1));
    dz = -dp./(RHOavg(i)*g);
    Ntheta = real(sqrt((g./RHOavg(i)).*(dptheta./dz)));
end
 
% solving for N = sqrt((-g/RHO)*(dRHO/dz))
 
for i = 2:n
    g = -9.8;
    RHOavg(i) = (RHO(i)+RHO(i-1))*0.5;
    dRHO = (RHO(i) - RHO(i-1));
    dz = -dp./(RHOavg(i)*g);
    NRHO = real(sqrt((g./RHOavg(i)).*(dRHO./dz)));
end
 
% solving for N = sqrt((-g/RHO)*(dRHOt/dz))
 
for i = 2:n
    g = -9.8;
    RHOavg(i) = (RHO(i)+RHO(i-1))*0.5;
    dRHOt = (RHOt(i) - RHOt(i-1));
    dz = -dp./(RHOavg(i)*g);
    NRHOt = real(sqrt((g./RHOavg(i)).*(dRHOt./dz)));
end
 
% graph of N vs Depth
 
figure(7);
plot(Ntheta,z);
xlabel('Buoyancy Oscillations (1/s)');
ylabel('Depth (m)');
hold on
plot(NRHO,z);
hold on 
plot(NRHOt,z);
legend({'RHO','RHOtheta','RHOt'},'Location','northeast');
set(gca,'YDir','reverse')

[Steven Lord]

What does "it is not working out" mean in this context?

• Do you receive warning and/or error messages? If so the full and exact text of those messages (all the text displayed in orange and/or red in the Command Window) may be useful in determining what's going on and how to avoid the warning and/or error.
• Does it do something different than what you expected? If so, what did it do and what did you expect it to do?
• Did MATLAB crash? If so please send the crash log file (with a description of what you were running or doing in MATLAB when the crash occured) to Technical Support using the Contact Support link on the Support section of the MathWorks website so we can investigate.

[Hannah Evans]

The graph shows up with no points on it, just an empty graph.

Answer 1

[Star Strider]

Subscript the variables in the loops:

Ntheta(i) = real(sqrt((g./RHOavg(i)).*(dptheta./dz)));

and similarly for the others in the plot calls.

I cannot run your code, so I am posting this as UNTESTED CODE. It should work.

[Hannah Evans]

It Worked! Thank you!!

[Star Strider]

My pleasure!