How to draw binary phase diagram ??

clear all;

clc;

%%

G1_s=300; % 1 element pure gibbs energy

G2_s=500; % 2 element pure gibbs energy

Om_s=5000; % 1,2 element mixing enthalpy

Om_l=3000; % 1,2 element mixing enthalpy

R=1.987; % gas constant

Tm_1=600; % element1 melting temparture

Tm_2=900; % element2 melting temparture

Hm_1=2000; % element1 Heat of fusion

Hm_2=1300; % element2 Heat of fusion

Tset=[0:700:1400]; %Temparture Setting

%%

for i =1:length(Tset)

T=Tset(i);

syms x a b G_l(x) G_s(x) dG_s(x) dG_l(x) Sola Solb

G1_l=G1_s+Hm_1/Tm_1*(Tm_1-T); % element1 liquid

G2_l=G2_s+Hm_2/Tm_2*(Tm_2-T); % element2 liquid

G_s(x) = G1_s*x+G2_s*(1-x)+Om_s*x*(1-x)+R*T*(x*log(x)+(1-x)*log(1-x)); % x-G function of Gibbs energy of solid

G_l(x) = G1_l*x+G2_l*(1-x)+Om_l*x*(1-x)+R*T*(x*log(x)+(1-x)*log(1-x)); % x-G function of Gibbs energy of liquid

dG_s(x)=diff(G_s(x),x);

dG_l(x)=diff(G_l(x),x);

denk1=(((G_s(b)-G_l(a))/(b-a))==(dG_s(a))); %eq 1

denk2=(dG_s(b)==dG_l(a)); %eq 2

[Sola,Solb] = vpasolve([denk1 denk2],[a b],[0.1;0.9]);

a=Sola;

b=Solb;

disp([T,a,b])

end

i think there is something wrong in my code and i want to ask you what is wrong

i want to plot T versus a graph and T versus b

Answer 1

[SaiDileep Kola]

Since the solution from vpasolve is complex in this case, try to plot as real vs complex values as described here.