function [E10 M0] = SI2con_fit(xdata, ydata, TR, M00) 

%global E1_weight;

    E1_weight = 1000;
    %par0 = [2e4 E1_weight*exp(-TR/500)];
    par0 = double([M00 E1_weight*exp(-TR/500)]);
    
%     options=optimset('TolX',1e-10,'TolFun',1e-20,'Display','iter','DiffMinChange',1e-5,...
%         'MaxFunEvals',1e10,...
%         'DerivativeCheck','off','GradObj','off','LargeScale','off',...
%         'MaxIter',50,'TolCon',1e-20);%,'PlotFcns',@optimplotfval,'PlotFcns',@optimplotx
    %[par,fval,exitflag,output,lambda,grad] = fmincon(@SI2con_precontrast,par0,[],[],[],[],[0 0],[1e10 E1_weight])
    %par = fminunc(@SI2con_precontrast,par0)

    %options = optimset('Jacobian','on');
    %par = lsqcurvefit(@SI2con_precontrast,par0, xdata, ydata); %, [0 0], [1e10 E1_weight], options)
    
    options=optimset('Display','off');
    myfun = @(par,xdata)SI2con_precontrast(par,xdata,E1_weight);
    par = lsqcurvefit(myfun,par0, xdata, ydata, [0 0], [1e10 E1_weight],options);
    %par = lsqcurvefit(@SI2con_precontrast,par0, xdata, ydata, [0 0], [1e10 E1_weight]);
    %par = par0;

    M0 = par(1);
    E10 = par(2)/E1_weight;
