Fsolve for multiple coupled-equation

Hi Guys, I have coupled-equation generated as T increased from 300 to 1300. Each 1 degree of T produce 1 coupled-equation I need to solve, so I have total 1000 coupled equation. Please help me to insert fsolve code into my loop. Thank you.

Note: I have found x0, written as Gmd. and here's the code:

 syms x
 %R,Regnault constant (gas constant) = 8.3144598 J/(K.mol)
 R=8.3144598
 %temperature
 T=300
 Z=1

 while T<=1300

    %Ga1, Gibbs standard for Mg alpha(hcp)
    if T>=298 & T<=923
        Ga1(Z)= (-8367.34)+(143.675547.*T)-(26.1849782.*T.*log(T))+(0.4858*(10.^-3).*(T.^2))-(1.393669*10.^-6.*T.^3)+   (78950.*T.^-1);
    else T>923 & T<=1300
        Ga1(Z)=(-14130.185)+(204.716215.*T)-(34.3088.*T.*log(T))+(1038.192.*10.^25.*T.^-9);
    end

    %Gb1, Gibbs standard for Li alpha(hcp)
    if T>=200 & T<=453
        Gb1(Z)=(-10737.817)+(219.637482.*T)-(38.940488.*T.*log(T))+(35.466931.*10.^-3.*T.^2)-(19.869816.*10.^-6.*T.^3)+(159994.*T.^-1);
    elseif T>453 & T<=500
        Gb1(Z)=(-559733.123)+(10549.879893.*T)-(1702.8886493.*T.*log(T))+(2258.3294448.*10.^-3.*T.^2)-(571.066077.*10.^-6.*T.^3)+(33885874.*T.^-1);
    else T>500 & T<=1300
        Gb1(Z)=(-9216.994)+(181.278285.*T)-(31.2283718.*T.*log(T))+(2.633221.*10.^-3.*T.^2)-(0.438058.*10.^-6.*T.^3)-(102387.*T.^-1);
    end

    %Ga2, Gibbs standard for Mg beta(bcc)
    if T>=298 & T<=923
        Ga2(Z)=(-5267.34)+(141.575547.*T)-(26.1849782.*T.*log(T))+(0.4858.*10.^-3.*T.^2)-(1.393669.*10.^-6.*T.^3)+(78950.*T.^-1);
    else T>923 & T<=1300
        Ga2(Z)=(-11030.185)+(202.616215.*T)-(34.3088.*T.*log(T))+(1038.192*10.^25.*T.^-9);
    end

    %Gb2, Gibbs standard for Li beta(bcc)
    if T>=200 & T<=453
        Gb2(Z)=(-10583.817)+(217.637482.*T)-(38.940488.*T.*log(T))+(35.466931.*10.^-3.*T.^2)-(19.869816.*10.^-6.*T.^3)+(159994.*T.^-1);
    elseif T>453 & T<=500
        Gb2(Z)=(-559579.123)+(10547.879893.*T)-(1702.8886493.*T.*log(T))+(2258.329444.*10.^-3.*T.^2)-(571.066077.*10.^-6.*T.^3)+(33885874.*T.^-1);
    else T>500 & T<=1300
        Gb2(Z)=(-9062.994)+(179.278285.*T)-(31.2283718.*T.*log(T))+(2.633221.*10.^-3.*T.^2)-(0.438058.*10^-6.*T.^3)-(102387.*T.^-1);
    end

    %Xb=mol fraction Li, Xa=mol fraction Mg, x defined as Xb so Xa=(1-x)

    Gm1= (1-x).*Ga1+x.*Gb1+R.*T.*((1-x).*log(1-x)+x.*log(x))+(1-x).*x.*((-6856)+4000.*(1-2.*x)+4000.*(1-2.*x).^2)
    Gm2= (1-x).*Ga2+x.*Gb2+R.*T.*((1-x).*log(1-x)+x.*log(x))+(1-x).*x.*((-18335)+8.49.*T+3481.*(1-2.*x)+(2658-0.114.*T).*(1-2.*x).^2)

    %assumption Gma=Gmb to find roots x0, intersection two phase alpha-beta
    %G1 is difference of Gibbs Li and Mg
    Gm=Gm2-Gm1;

    %roots of Gm, Gmd=x0
    %Gmx first derrivative, Gmy second derivative

    for i=1:length(Gm)
        Gma=Gm(i);
        Gma==0
        Gmb=solve(Gma,x);
        Gmd(i)=Gmb(2)
    end

    %fsolve to find x1 and x2

    % loop for next T value
    T=T+1
    Z=Z+1
 end