Calling a function iteratively within a for loop with various size outputs

Hi this could be a long one, thanks for your consdieration. I have a "for" loop which iterates over a growing differential equation (the inputs of thuis may or may not be known a-priori). I then iterate over each extension of the differential equation and solve it. But I write the differential equation to a .m file which I then call as a function and solve. The output of this differential equation is a matrix N rows and M columns, The number of rows (1501) shouldn't change between each iteration however the number of columns M will. In this example it should grow in steps of 3, 6, 9, 12, 15. This could change later though so I wouldn't like them constrained. For some reason the code goes through the first iteration then at the line:

A=CME_solver(maxBeta);

The code breaks.

To explain the script I first decide what order of extension with the for loop. I then provide some logic where I select certain combinations of a vector.

I then add and subtract these elements of a vector to generate a matrix of sums and differences of each combination based on where each element came from in the original vector.

By comparing the numbers in the matrix with the numbers in the initial vector I can extract all the numbers (and locations) from the matrix for all numbes which existed in the initial vector and thus match what sum and difference combinations produced them.

I have to then convert this to a differential equation which is written as a string in two .m files with the names CME_solver.m and CME_solver_no_coeff.m. These Are standard ode45's essentially written programaticaly and then placed into the same folder as the script.

I then want to get the result of this differential equation and store each result for each iteration in the cell array I and I_no_coeff which i would finally pass back from this function to a main script once the for loop is finished.

Please see code below (uncomment the first inputs here to run as a stand alone script if necesary)

function [I,I_no_coeff,f_CME] = CMEj_solver_writer(param,N,maxBeta,Ap0,As0,up_folder)

% up_folder=/home/[your_folder_here]/..;
% Ap0=0.9502;
% As0=0.0801;
% param(5)=12;
% param(7)=7.2;
% N=1500;
% maxBeta=0.433;

CME_extension = 5;

for harmonic=1:CME_extension
    
    
    wp = 2*pi*param(5)*10^9;
    ws = 2*pi*param(7)*10^9;


    wi = wp-ws;

%    wdel = abs(ws-wi)

    w = [wp,ws,wi];

    w=sort(w)
    wc=harmonic*wp;


   for k=1:harmonic^2
        for j = 1:length(w)
            for i=1:length(w);
                a = w(j) + w(i);
                b = w(j) - w(i);
                check1 = ismembertol(a,w,0.001);
                check2 = ismembertol(b,w,0.001);
                if mod(w(j),wp)==0 || mod(w(i),wp)==0 && b~=0
                    if check1==0 && a<wc+1e6
                       w = [w,a];
                      
                    end
                    if check2==0 && b>0
                       w = [b,w];
                      
                    end
                    
                end
               
            end
        end
   end
    

    %toDelete=w>wc;
    %w(toDelete)=[];
    w=sort(w);
    f=w./(2*pi)
    fa=a./(2*pi)

    for j=1:length(w)
        for i=1:length(w)
            sfg = w(j) + w(i);
            check = ismembertol(sfg,w,0.01);
            if check==1;
                sfg_matrix(j,i) = sfg;
            else
                sfg_matrix(j,i) = 0;
            end
            if i>j
                sfg_matrix(j,i) = 0;
            
            end
            dfg = w(j) - w(i);
            check = ismembertol(dfg,w,0.01);
            if check==1;
                dfg_matrix(j,i) = dfg;
            else
                dfg_matrix(j,i) = 0;
            end
            
        end
    end


    dfg_trans = dfg_matrix';
    mixing_matrix = sfg_matrix + dfg_trans;

    for i = 1:length(w)
        w_string{i} = ['w',num2str(i)];
        k_string{i} = ['k',num2str(i)];
        amp_string{i} = ['A(',num2str(i),')'];
    end

    for j=1:length(w)
        for i=1:length(w)

            %%%% mixing_matrix(j,i) for j<i is dfg ie top right half
            %%%% mixing_matrix(j,i) for j<i is sfg ie bottom left half
            %%%% mixing_matrix(j,i) for j=i is SHG of ji = jj = ii is NW-SE
            %%%% diagonal
            if ismembertol(mixing_matrix(j,i),w,0.001)==1
                if j<i
                    [check,place] = ismembertol(mixing_matrix(j,i),w);
                    dA_string{place,i} = ['+ (',k_string{i},'-',k_string{j},')*',k_string{i},'*',k_string{j},'*',amp_string{i},'*conj(',amp_string{j},')*exp(1i*(',k_string{i},'-',k_string{j},')*x) '];
                end

                if j>i
                    [check,place] = ismembertol(mixing_matrix(j,i),w,0.001);
                    dA_string{place,i} = ['- (',k_string{i},'+',k_string{j},')*',k_string{i},'*',k_string{j},'*',amp_string{i},'*',amp_string{j},'.*exp(1i*(',k_string{i},'+',k_string{j},')*x) '];
                end

                if j==i 
                    if ismembertol(mixing_matrix(j,i),w,0.001)==1
                       [check,place] = ismembertol(mixing_matrix(j,i),w) ;
                       dA_string{i+j,i+j} = ['- 0.5*(',k_string{i},'+',k_string{j},')*',k_string{i},'*',k_string{j},'.*',amp_string{i},'*',amp_string{j},'.*exp(1i*(',k_string{i},'+',k_string{j},')*x) '];
                    end
                end
            end
        end
    end



    for j=1:length(w)
        for i=1:length(w)

            %%%% mixing_matrix(j,i) for j<i is dfg ie top right half
            %%%% mixing_matrix(j,i) for j<i is sfg ie bottom left half
            %%%% mixing_matrix(j,i) for j=i is SHG of ji = jj = ii is NW-SE
            %%%% diagonal
            if ismembertol(mixing_matrix(j,i),w,0.001)==1
                if j<i
                    [check,place] = ismembertol(mixing_matrix(j,i),w,0.001);
                    dA_string_no_coeff{place,i} = ['+ ',k_string{i},'*',k_string{j},'*',amp_string{i},'*conj(',amp_string{j},')*exp(1i*(',k_string{i},'-',k_string{j},')*x) '];
                end

                if j>i
                    [check,place] = ismembertol(mixing_matrix(j,i),w,0.0001);
                    dA_string_no_coeff{place,i} = ['- ',k_string{i},'*',k_string{j},'*',amp_string{i},'*',amp_string{j},'.*exp(1i*(',k_string{i},'+',k_string{j},')*x) '];
                end

                if j==i 
                    if ismembertol(mixing_matrix(j,i),w,0.0001)==1
                       [check,place] = ismembertol(mixing_matrix(j,i),w,0.001);
                       dA_string_no_coeff{i+j,i+j} = ['- 0.5*',k_string{i},'*',k_string{j},'.*',amp_string{i},'*',amp_string{j},'.*exp(1i*(',k_string{i},'+',k_string{j},')*x) '];
                    end
                end
            end
        end
    end




    for i=1:length(dA_string)
       dA_string_final{i} = ['(maxBeta/(2*',k_string{i},'))*(',dA_string{i,:},')*exp(-1i*',k_string{i},'*x);'];

    end

    for i=1:length(dA_string)
       dA_string_no_coeff_final{i} = ['(maxBeta/(2))*(',dA_string_no_coeff{i,:},')*exp(-1i*',k_string{i},'*x);'];

    end



    fid = fopen([up_folder,'/CME_solver.m'],'w');

    fprintf(fid, 'function [A] = CME_solver(maxBeta) \r\n\r\n');

    for i=1:length(w)
       k=acos(1-((w.^2).*param(3).*1e-12.*param(2).*1e-15)./(2.*(1-(w.^2).*param(3).*1e-12.*param(8).*1e-15)));
       fprintf(fid, '%s = %0.4f\r\n',k_string{i},k(i));
    end

    fprintf(fid, '\r\n x=[0:%i] \r\n\r\n',N-1)

    for i=1:length(dA_string_final)
        if i==1
            fprintf(fid, 'dA5 = @(x,A)[%s \r\n',dA_string_final{i});
        elseif i==length(dA_string_final)
            fprintf(fid, '%s] \r\n',dA_string_final{i});
        else
            fprintf(fid, '%s \r\n',dA_string_final{i});
        end
    end

    fprintf(fid, '\r\n [x,A] = ode45(dA5, x, [')
    for i =1:length(w)
        if ismembertol(w(i),ws,0.001)==1
            fprintf(fid, '%s ',As0)
        elseif ismembertol(w(i),wp,0.001)==1
            fprintf(fid, '%s ',Ap0)
        else 
            fprintf(fid, '0 ')
        end
    end
    fprintf(fid, '])')
    fclose(fid)


    fid = fopen([up_folder,'/CME_solver_no_coeff.m'],'w')

    fprintf(fid, 'function [A] = CME_solver_no_coeff(maxBeta) \r\n\r\n')

    for i=1:length(w)
       k=acos(1-((w.^2).*param(3).*1e-12.*param(2).*1e-15)./(2.*(1-(w.^2).*param(3).*1e-12.*param(8).*1e-15)));
       fprintf(fid, '%s = %0.4f\r\n',k_string{i},k(i)) 
    end

    fprintf(fid, '\r\n x=[0:%i] \r\n\r\n',N-1)

    for i=1:length(dA_string_final)
        if i==1
            fprintf(fid, 'dA5 = @(x,A)[%s \r\n',dA_string_no_coeff_final{i})
        elseif i==length(dA_string_final)
            fprintf(fid, '%s] \r\n',dA_string_no_coeff_final{i})
        else
            fprintf(fid, '%s \r\n',dA_string_no_coeff_final{i})
        end
    end

    fprintf(fid, '\r\n [x,A] = ode45(dA5, x, [')
    for i =1:length(w)
        if ismembertol(w(i),ws,0.001)==1
            fprintf(fid, '%s ',As0)
        elseif ismembertol(w(i),wp,0.001)==1
            fprintf(fid, '%s ',Ap0)
        else 
            fprintf(fid, '0 ')
        end
    end
    fprintf(fid, '])')
    fclose(fid)

    %clear dA_string dA_string_final dA_string_no_coeff_final dA_string_no_coeff
    Lg = param(3)*10^-12;
    Ic = param(10)*10^-6;
    Cg = param(2)*10^-15;
    Cj = param(8)*10^-15;
    phi0 = 2.067833848*10^-15;
    Bl = 2*pi*Lg*Ic/phi0;


    A=CME_solver(maxBeta);

    V=(Lg*Ic/Bl).*w.*abs(A);

    Z=sqrt(Lg/Cg);

    I{harmonic} = V/Z;

    A=CME_solver_no_coeff(maxBeta);

    V=(Lg*Ic/Bl).*w.*abs(A);

    Z=sqrt(Lg/Cg);

    I_no_coeff{harmonic} = V/Z;

    f_CME{harmonic}=w./(2*pi);
end

Thankyou for your consideration, it is much appreciated.

[Sindar]

What is the actual error given?

Without looking too closely at your code, is it possible that the new version of the function has not been processed by Matlab, so it is using the old code where the size doesn't fit? This might be difficult to work around; is it possible you could generalize the function in such a way that it can take inputs to "update" rather than programmatically writing code?

[Thomas Dixon]

Hi Sindar. Sorry for late response. I actually got this to work, but now I want to do exactly as you suggest. I dont want to have text files of pseudocode written and then read, as this seems fundementally stupid. (the worst kind of stupid). I want to somehow create the dA vector programaticaly based on the input A vector. I really am struggling with this at the moment. In the meantime I have created an app with default values which will calculate the mixing matrix and hopefully give insight into my problems.

You have to click 'Calculate' next to BetaL ( you may wish to click on the plot button as well to see my wonderful plotting)

Then go to the next tab ( I recommend that you choose MLine from the Cut-off determination tab and flick between All tones and Pump mediated for the Mixing mehcnaisms tab. The red light on Calculating tones means the script is running. I suggest for simplicity sticking with Pump mediated as it will be equal to or les than the number of entries in all tones.

Finally hit the Calculate Mixing Matrix buton and a table wil appear which has values corresponding with the values in the frequencies table (top right).

Whew.....

Now when you look at the mixing matrix table called mixingmatrix(i,j) then there is an obvious rule of correspondence between the numbers in Frequency=F(n) and MixingMatrix(i,j) the rule is:

for i=1:length(F)
    for j=1:length(F)
        
        
        if i>j
            mixingmatrix(i,j)=F(i)-F(j)
        end
        
        if i<=j
            mixingmatrix(i,j)=F(i)+F(j)
        end

    end
    
end

These mixing conditions and how theye are created define the relationship between A(n) and dA(n) which is always the same length of F(n). So with this rule and indexing I can work out the functional form of dA with regard to A. The problem I have is programtically giving it to matlab. I just don't know how to do it.

Thanks again for any help.