Need help modifying program to run faster

Hi, I've developed the following program and sub-program that essentially helps me with fitting data to a specified logistic function and analyzing the returned regression coefficients and their corresponding standard errors. The problem is that the program takes a very long time to run when larger data sets are used. I've tried to make the program faster in everyway I know how (i.e. matrix pre-allocation, function handles, etc.) but with my limited knowledge of Matlab I don't know what else I can do to speed things up. Are there any modifications that I can make to the main program or the sub-program that would help speed it up? Also, since the program spends most of its time running through the sub-program 'nlinfit' in the for-loop, does anyone have any ideas of how I could fit the data to the logistic function on my own istead of using nlinfit? - Any help would be greatly appreciated, thanks.

if true
  % code
end
%------------------------------------------------------------------------
% Gaussian_Beam_C_Coefficient_analysis.m
%
% PROGRAM DESCRIPTION:
%   1.) Reads in a set of data points from a user-specified Excel file
%   2.) Fits the data to a nonlinear logistic model and computes the
%        the regression coefficients by using a nonlinear least-squares 
%        regression method in the function nlinfit.m
%   3.) Computes the standard error for the regression coefficients
%        by using the function StdErrFunc.m, a modified version of
%        the function nlparci.m
%   4.) Uses the function nchoosek.m to compute all possible
%        combinations with four points from the data set, fits each 
%        combination to the same nonlinear logistic model, and computes
%        the regression coefficients and their corresponding standard
%        errors for each combination.
%   5.) Writes the specified data to an Excel file
%
% Input:
%   Data     -   Matrix of data points from Excel file
%   [cc,rc]  -   Column count and row count of data matrix
%   n_o_p    -   Number of points in each combination
%   n_o_c    -   Number of combinations possible in the data set
%   xD       -   Column of x-values from matrix 'Data'
%   yD       -   Column of y-values from matrix 'Data'
%   x1       -   Row of x-values
%   IRCP     -   Initial regression coefficient predictions
%   CI       -   Confidence interval
%
% Computations:
%   LM1      -   Logistic model
%                 Form:                    1/(RC(1)+RC(2)*e^(-RC(3)*x))
%                 Regression coefficients: RC(1) = A
%                                          RC(2) = B
%                                          RC(3) = C
%   [RC,R,J] -   Parameters of logistic model
%                 RC    = Regression coefficients
%                 R     = Residuals
%                 J     = Jacobian of specified logistic model
%   se1      -   Standard error values of regression coefficients
%   iA       -   Initial value of A-coefficient
%   iAse     -   Initial standard error of A-coefficient
%   iB       -   Initial value of B-coefficient
%   iBse     -   Initial standard error of B-coefficient
%   iC       -   Initial value of C-coefficient
%   iCse     -   Initial standard error of C-coefficient
%   xD_combs -   All possible combinations of x-values from the data set
%   yD_combs -   All possible combinations of y-values from the data set
%   n        -   number of loops
%   M        -   Matrix that contains the x- and y- values of the nth
%                 combination
%   xM       -   Column of x-values from matrix 'M'
%   yM       -   Column of y-values from matrix 'M'
%   x2       -   Row of x-values from matrix 'M'
%   se2      -   Standard error values of regression coefficients of
%                 nth combination
%   A        -   A-coefficient of nth combination
%   Ase      -   Standard error of A-coefficient of nth combination
%   B        -   B-coefficient of nth combination
%   Bse      -   Standard error of B-coefficient of nth combination
%   C        -   C-coefficient of nth combination
%   Cse      -   Standard error of C-coefficient of nth combination
%   MM       -   Matrix containing recomputed parameters of C-coefficient
%
% Output:
%   Writes the data contained in matrix MM to a specified Excel file.
%
% Written by: Christopher Thomas Collinsworth
% Written on: September 15, 2013
%------------------------------------------------------------------------
%%Clear variables and command window
clear all
close all
clc
%%Read in Excel data----------------------------------------------------------------
Data = xlsread('GBD.xlsx');
%%Input-----------------------------------------------------------------------------
[rc,cc]   = size(Data);
n_o_p     = 4;                                          % User defined
n_o_c     = nchoosek(rc,n_o_p);
xD        = Data(:,1);
yD        = Data(:,2);
x1        = xD';
IRCP      = [0.019 221.8 3.8]';
CI        = 0.05;
%%Calculate initial parameters------------------------------------------------------
LM1       = @(RC,xD)(1./(RC(1)+RC(2)*exp(-RC(3)*x1)))'; % function handle
[RC,R,J]  = nlinfit(xD,yD,LM1,IRCP);                    % fits data to model
se1       = StdErrFunc(RC,R,'jacobian',J,'alpha',CI);   % computes standard error
iA        = RC(1,1);
iAse      = se1(1,1);
iB        = RC(2,1);
iBse      = se1(2,1);
iC        = RC(3,1);
iCse      = se1(3,1);
%%For-loop--------------------------------------------------------------------------
xD_combs  = nchoosek(xD(1:rc),n_o_p)';
yD_combs  = nchoosek(yD(1:rc),n_o_p)';
MM        = zeros(n_o_c,7);                             % Matrix pre-allocation
tic
for n     = 1:1:n_o_c
M        = [xD_combs(1:n_o_p,n) yD_combs(1:n_o_p,n)];
xM       = M(:,1);
yM       = M(:,2);
x2       = xM';
LM2      = @(RC,xM)(1./(RC(1)+RC(2)*exp(-RC(3)*x2)))'; % function handle
[RC]     = nlinfit(xM,yM,LM2,IRCP,CI);                 % fits data to model
se2      = StdErrFunc(RC,R,'jacobian',J,'alpha',CI);   % computes standard error
A        = RC(1,1);
Ase      = se2(1,1);
B        = RC(2,1);
Bse      = se2(2,1);
C        = RC(3,1);
Cse      = se2(3,1);
MM(n,:)  = [n A Ase B Bse C Cse];
end %for
%%Write data to Excel file----------------------------------------------------------
xlswrite('GBDTESTOUTPUT.xlsx',MM)
toc

if true
  % code
end
%------------------------------------------------------------------------------------
%==== StdErrFunc ====
function se = StdErrFunc(beta,resid,varargin)
J     = [];
alpha = 0.05;
if nargin>=3 && ischar(varargin{1})
  % Calling sequence with named arguments
  okargs =   {'jacobian' 'covariance' 'alpha'};
  defaults = {[]         []           0.05};
  [J, Sigma, alpha] = internal.stats.parseArgs(okargs,defaults,varargin{:});
end
% Remove missing values.
resid = resid(:);
n     = length(resid);
p     = numel(beta);
v     = n-p;

 [~,R] = qr(J,0);
 Rinv = R\eye(size(R));
 diag_info = sum(Rinv.*Rinv,2);
 rmse = norm(resid) / sqrt(v);
 se = sqrt(diag_info) * rmse;
end

[Jan]

As usual: No, the action performed by clear all are not correctly described by "clear variables". This removes all loaded M- and MEX-files from the memory also and reloading can waste seconds. In addition all breakpoints of the debugger are removed and everything, which impedes debugging in general is a really bad idea. Perhaps you mean clear variables, which is less brute.