function hd=ExoterDE(M0,sigmaMun,R0,sigmaRun)
%%Constant:
G=6.67428e-11;
Mearth=5.9736e24;
Rearth=6372797;

step=0.005;

%find how much mass uncertainty is equivalent to radius uncertainty in the
%sense of composition
x=[0.5,1,2,4,8,16,20];
y=[3.3,3.4,3.6,3.9,4.1,4.4,4.5];
ratio=interp1(x,y,M0,'pchip','extrap');
correction=sqrt((sigmaMun/M0)^2+(ratio*(sigmaRun/R0))^2)/(abs(sigmaMun/M0)+ratio*abs(sigmaRun/R0));
sigmaM=correction*sigmaMun;
sigmaR=correction*sigmaRun;

%%create the picture
figure;
colormap;
terplot(11);
terlabel('\beta (silicate)','\alpha (iron)','\gamma (water)');
title(['Mass:',num2str(M0),'M_\oplus',' \sigma_M:',num2str(sigmaMun),'M_\oplus',...
       '   Radius:',num2str(R0),'R_\oplus',' \sigma_R:', num2str(sigmaRun),'R_\oplus'],'FontSize',14,'FontWeight','bold');
hold on;

%%Data of EOS:
A=load('iron.txt');
B=load('silicate.txt');
C=load('water.txt');

%%Mass, Radius and their uncertainties we know:

for v=-3:3
M=M0*(1+v*sigmaM/M0)*Mearth;
R=R0*(1-v*sigmaR/R0)*Rearth;


%%Initial Conditions:
m0=M;
P0=0;

%%control parameters:
i=0;
j=0;
k=0;
w(1:7)=0;
a=0;


while ((a>=0)&(a<=1))&(j*k<0.5)
    mid=(1-a)/2;
    b=mid;
    l=mid/2;
    u1=0;     %%control parameter
    u2=0;     %%comtrol parameter
    while (l>5e-4)
        %%Solve the ODEs:
        options = odeset('RelTol',1e-4,'AbsTol',[m0*1e-4 3.6e7]);
       [d,y]=ode45('Exoeqnew',[0 0.9999*R],[m0,P0],options,G,a,b,m0,R,A,B,C);
        sign=y(length(y),1)/M;
        %%Improve to Newton's method
        if sign>0
            b=b+l;u1=1;
        else
            b=b-l;u2=1;
        end
    l=l/2;
    clear y;
    end
    if u1*u2>0.5
        i=i+1;iron(v+4,i)=a;silicate(v+4,i)=b;j=1;k=0;
    else
        k=1;
    end
    a=a+step;
end

if i>3
    p1=polyfit(iron(v+4,1:i),silicate(v+4,1:i),3);
    silicatenew(v+4,1:i)=polyval(p1,iron(v+4,1:i));
    hd(v+4)=ternplot(iron(v+4,1:i),silicatenew(v+4,1:i),round(abs(v)+1));
    hold on;
elseif i>0.5
    hd(v+4)=ternplot(iron(v+4,1:i),silicate(v+4,1:i),round(abs(v)+1));
    hold on;
else
    hd(v+4)=0;w(v+4)=1;
    ['curve for Mass=',num2str(M0*(1+v*sigmaM)),'M_\oplus',' Radius=',num2str(R0*(1-v*sigmaR)),'R_\oplus',...
        ' does not exist.'];
end
end
hold off;