How to deal properly with a fitting to a non-linear curve nested in a for loop

Hello everyone,

I am dealing with the following situation. I have the following arrays in my workspace: data_set_elastic (30001x4x2431 elements), physical_time_elastic (2431x1 elements), and spatial_grid (30001x1 elements). Each index of the third dimension of data_set_elastic defines a plane 30001x4 where the first column represents spatial position (which is exactly the same that it is in spatial_grid), and the rest of them magnetization components, mx, my and mz, correspondingly. For each index of the third dimension I can know which is the related time instant to it, because I have an array physical_time_elastic (2431x1 elements). So, at the end, I have encoded in data_set_elastic which are the values for mx, my and mz for each time and spatial position. Now, for each instant of time (that it is, for each plane 30001x4 defined by each element of the third dimension of the array data_set_elastic) I want to fit the data in the second element of the second dimension of data_set_elastic. That it is, I will fit for each instant of time, along all the spatial dimension, the magnetization profile in the x-th direction (second element of the second dimension of the array data_set_elastic). To simplify my problem, I have created a two-dimensional array for mx only:

mx_elastic=zeros(length(spatial_grid),length(physical_time_elastic));
for i=1:length(physical_time_elastic)
    for j=1:length(spatial_grid)
        mx_elastic(j,i)=data_set_elastic(j,2,i);
    end
end

The function to which I want to fit my data mx_elastic (30001x2431 elements) is the following fitting_DWs_profile_elastic (the rest of the information is the creation of the arrays that will encode the fitting parameters A1_elastic, A2_elastic, q1_elastic, q2_elastic, Delta1_elastic, and Delta2_elastic, and a way to find a proper initial value for each of these fitting parameters):

A1_elastic=zeros(length(physical_time_elastic),1);
A2_elastic=zeros(length(physical_time_elastic),1);
q1_elastic=zeros(length(physical_time_elastic),1);
q2_elastic=zeros(length(physical_time_elastic),1);
Delta1_elastic=zeros(length(physical_time_elastic),1);
Delta2_elastic=zeros(length(physical_time_elastic),1);
[median_value_mx_elastic,index_median_value_mx_elastic]=...
    min(abs(data_set_elastic(:,1,1)-median(data_set_elastic(:,1,1))));
[max_mx_value_elastic_left_part,index_max_mx_value_elastic_left_part]=...
    max(data_set_elastic((1:(index_median_value_mx_elastic-1)),2,1));
zci_left_part_elastic=@(v_left_part_elastic) find(v_left_part_elastic(:)...
    .*circshift(v_left_part_elastic(:),[-1 0])<=0); 
mx_temp_left_part_elastic=data_set_elastic((1:(index_median_value_mx_elastic-1)),2,1)...
    -0.5*max(data_set_elastic((1:(index_median_value_mx_elastic-1)),2,1));
mx_dx_left_part_elastic=zci_left_part_elastic(mx_temp_left_part_elastic);
Delta1_elastic(1)=diff(spatial_grid(mx_dx_left_part_elastic));
[min_mx_value_elastic_right_part,index_min_mx_value_elastic_right_part]=...
    min(data_set_elastic(((index_median_value_mx_elastic+1):end),2,1));
A1_elastic(1)=max_mx_value_elastic_left_part;
A2_elastic(1)=min_mx_value_elastic_right_part;
q1_elastic(1)=data_set_elastic(0+index_max_mx_value_elastic_left_part,1,1);
q2_elastic(1)=data_set_elastic(index_median_value_mx_elastic+...
    index_min_mx_value_elastic_right_part,1,1);
Delta1_elastic(1)=fwhm(spatial_grid,data_set_elastic(1:...
    (index_median_value_mx_elastic-1),2,1));
Delta2_elastic(1)=fwhm(spatial_grid,data_set_elastic(index_median_value_mx_elastic+1:...
    end,2,1));
fitting_DWs_profile_elastic=@(A1_elastic,q1_elastic,Delta1_elastic,A2_elastic,...
    q2_elastic,Delta2_elastic)(A1_elastic.*sech((spatial_grid-q1_elastic)./Delta1_elastic)+...
    abs(A2_elastic).*sech((spatial_grid-q2_elastic)./Delta2_elastic));

Once I have this, I want to implement a for loop to do the fitting of my data encoded in mx_elastic to the fitting_DWs_profile_elastic:

for i=2:length(physical_time_elastic)
    fitted_DW_elastic=fit(fitting_DWs_profile_elastic,spatial_grid,mx_elastic(:,i),...
        [A1_elastic(i-1),A2_elastic(i-1),q1_elastic(i-1),q2_elastic(i-1),Delta1_elastic(i-1),...
        Delta2_elastic(i-1)]);
    A1_elastic(i)=fitted_DW_elastic.A1_elastic;
    A2_elastic(i)=fitted_DW_elastic.A2_elastic;
    q1_elastic(i)=fitted_DW_elastic.q1_elastic;
    q2_elastic(i)=fitted_DW_elastic.q2_elastic;
    Delta1_elastic(i)=fitted_DW_elastic.Delta1_elastic;
    Delta2_elastic(i)=fitted_DW_elastic.Delta2_elastic;
end

This does not work, and MatLab gives me the following message:

Undefined function 'fit' for input arguments of type 'function_handle'.

Moreover, I want to impose an additional condition inside the for loop. When q1_elastic(i)=q2_elastic(i), stop the fitting, and even do not give me the values of the fitting for that i index. Moreover, this must imply also that the dimension of the A1_elastic, A2_elastic, q1_elastic, q2_elastic, Delta1_elastic, and Delta2_elastic must be length(i for which q1_elastic=q2_elastic). However, as you see, I have imposed at the beginning that the length of the array of these fitting parameters it is equal to the number of elements evaluated in the for loop. An attempt without knowing too much how this will work:

for i=2:length(physical_time_elastic)
    fitted_DW_elastic=fit(spatial_grid,mx_elastic(:,i),fitting_DWs_profile_elastic,...
        [A1_elastic(i-1),A2_elastic(i-1),q1_elastic(i-1),q2_elastic(i-1),Delta1_elastic(i-1),...
        Delta2_elastic(i-1)]);
    A1_elastic(i)=fitted_DW_elastic.A1_elastic;
    A2_elastic(i)=fitted_DW_elastic.A2_elastic;
    q1_elastic(i)=fitted_DW_elastic.q1_elastic;
    q2_elastic(i)=fitted_DW_elastic.q2_elastic;
    Delta1_elastic(i)=fitted_DW_elastic.Delta1_elastic;
    Delta2_elastic(i)=fitted_DW_elastic.Delta2_elastic;
    if q1_elastic(i)=q2_elastic(i)
         break
    else continue
    end
end

Any suggestion?

If I check which -all fit, I found the following:

>> which -all fit
C:\Program Files\MATLAB\R2019b\toolbox\stats\stats\@gmdistribution\fit.m  % gmdistribution method