Hello,
I did some modifications to your code and I think I got it most of the way there but it is hard to finish because I am not quite sure how to structure the data for subplot 2 since I do not know how the variables relate to each other.
One thing I did was vectorize the results of your calculations and run the plots only once, that sped things up and seemed to clean up the second plot. Keep in mind that when you have a plot command with "hold on" in a loop, not only does MATLAB call the plot function on every loop cycle, which can cause high processing times, it can also be very difficult to troubleshoot when you don't get the results you expect.
I think I got you most of the way there but I was not sure how to structure the x-axis data for the second subplot so I just created a generic index vector so I could look at the shape of the plot and see if it was still doing crazy things.
Hope this at least helps:
%Written by Adam Kelly
clear all;
clc;
%constants
r=1.5e+10;
rs=6.261e+8; %Radius of Star
ro=1.49e11; %Sun orbital radius (1 AU)
rRD=1.25e8 ; %Radius of Red Dwarf
stefconst=5.67e-8; %stefan-boltzmann constant
T=5200.2; %Temp of Sun
TRD=3773.15; %Temp of Red Dwarf
al=0.3; %albedo
a=3; %Red Dwarf's semi major axis
G=6.67408e-11;
MS=2e+30;
MP=5.972e+24;
MRD=1.2e+30;
%luminosity of star
Ls=4*pi*rs^2*stefconst*T^4;
Lrd=4*pi*rRD^2*stefconst*TRD^4;
fluxS=(Ls./(4*pi*ro^2));
%Red Dwarf's Period
%Using Kepler's 3rd law
P=a.^1.5; %Gives answer in years
PRD=P*(365)
vTP=[];
roRD = 3.75:-0.01:0.25;
for i = 1:length(roRD)
%Flux
fluxRD=(Lrd/(4*pi*(roRD(i)*1.496e11)^2)); %flux of red dwarf
fluxAverage=(fluxS+fluxRD);
%Temp of planet
Tplanet=((fluxAverage*(1-al))/(4*stefconst))^0.25;
TP(i)= Tplanet-273.15;
vTP(i) = TP(i);
end
subplot(2,1,1);
plot(roRD,TP,'.r')
xlim([0 4])
ylim([-58 -44])
xlabel('Orbit(AU)')
ylabel('Temperature of Planet(C)')
plotVectorX = 1:length(vTP);
subplot(2,1,2);
plot(plotVectorX,vTP,'.b');
xlim([0 1896])
ylim ([-58 -45])
