Are there any prebuilt operations to optimize the fit of a simulation to experimental data?

Function CVsims()

// VARIABLES I WANT TO OPTIMIZE

Variable C = 1E-6;      // moles/cm^3, Initial concentration of R.
Variable D = 8E-7;      // cm^2/s, Diffusion coefficient
Variable k0 = 5E-2;     // cm/s, electrochemical rate constant.

// OTHER VARIABLES THAT ARE EXPERIMENT DEPENDENT BUT FIXED

Variable Ehigh = 0.6;   // V, high potential of scan
Variable Elow = -0.1;   // V, low potential of scan
Variable E0 = 0.25;     // V, formal reduction potential of the compound
Variable n = 1;         // unitless, number of electrons transferred
Variable alpha = 0.1;   // unitless, alpha
Variable kc = 0;            // 1/s, rate constant for degredation of species R
Variable A = pi*(0.3/2)^2; // cm^2, Electrode area

// PHYSICAL CONSTANTS
Variable T = 298.15;        // Temp

Variable F = 96485;     // C/mol, Faraday's constant
Variable R = 8.3145;        // J/(mol*K) Ideal gas constant
Variable f_ = F/(R*T);  // 1/V, normalized Faraday's constant to room temp.

// SIMULATION VARIABLES
Variable L = 2000;      // Number of time points in the scans. Larger numbers give smoother CVs, but take longer.
Variable DM = 0.45;     // Model diffusion coefficient (keep below 0.5)
Variable j = ceil(4.2*sqrt(L))+5; // Number of distance boxes from electrode surface.

// Specify the scan rates to simulate (in V/s)
// Note that if you change N, it will alter the number of current waves.
// You'll need to re-adjust the plot.
Make /O /N=(6) scan;    
scan = {0.001, 0.005, 0.01, 0.050, 0.1, 1};

// This makes a 1-d array of the following variables as a function of scan rate.
// tk us the total time of the CV (based on the potential window and scan rate)
Make /O /N=(numpnts(scan)) tk = 2*(Ehigh-Elow)/scan;    // units = s
// Dt is the time length of each time step.
Make /O /N=(numpnts(scan)) Dt = tk/L; // units = s
// Dx is the width of each distance step from the electrode surface.
Make /O /N=(numpnts(scan)) Dx = sqrt(D*Dt/DM); // units = cm
// ktk is the dimensionless kinetic parameter
Make /O /N=(numpnts(scan)) ktk = kc*tk;
// km is the normalized diemnsionless kinetic parameter
Make /O /N=(numpnts(scan)) km = ktk/L;

Variable h1,h2; // These variables are used for the following for lop.

Make /O /N=(L+1,numpnts(scan)) k; // k is the time index (time point)
Make /O /N=(L+1,numpnts(scan)) t_;// t_ is the simulation time in seconds.

for (h1=0; h1<(L+1); h1++) // scanning through all time points.
    for (h2=0; h2<(numpnts(scan)); h2++) // scanning through each scan rate
   
        k[h1][h2] = h1; // this sets the time point number.
        t_[h1][h2] = Dt[h2]*k[h1][h2]; // this sets the actual amount of time at each time point number as a function of scan rate.
       
    endfor
endfor


Make /O /N=(L/2) E1 = Elow + 2*(Ehigh-Elow)*x/L; // This is the potential of the first scan direction at each time step.
Make /O /N=(L/2+1) E2 = Ehigh - 2*(Ehigh-Elow)*x/L; // This is the potential of the second scan direction at each time point.
Concatenate /NP /O {E1,E2},Escan; // This combines the potentials of the forward and reverse scans.
Make /O /N=(L+1) Enorm = Escan*f_; // This normalizes the calculated potential to the unitless value.
Make /O /N=(L+1) kf = k0*exp(-alpha*n*(Enorm-E0*f_)); // The forward rate constant as a function of potential.
Make /O /N=(L+1) kr = k0*exp((1-alpha)*n*(Enorm-E0*f_)); // The reverse rate constant as a function of potential.

Make /O /N=(L+1,j,numpnts(scan)) Ox = 0; // The initial concentration of O in solution normalized to 1.  
Make /O /N=(L+1,j,numpnts(scan)) Red = 1; // The initial concentration of R in solution normalized to 1.
Make /O /N=(L+1,numpnts(scan)) JO = 0; // The flux matrix.
Make /O /N=(l+1) Current;                    // A temporary holder to calculate the current from the flux.

Variable i1,i2,iv;

for (iv=0; iv<numpnts(scan); iv++) // for-loop for different scan rates.

    for (i1=0; i1<l; i1++) // for-loop for time steps.
   
        for (i2=2; i2<j-1; i2++) // for-loop for distance steps.
       
            // These represent diffusion from cell to cell.
            Ox[i1+1][i2][iv] = Ox[i1][i2][iv] + DM*(Ox[i1][i2+1][iv]+Ox[i1][i2-1][iv]-2*Ox[i1][i2][iv]);
            Red[i1+1][i2][iv] = Red[i1][i2][iv] + DM*(Red[i1][i2+1][iv]+Red[i1][i2-1][iv]-2*Red[i1][i2][iv]) - km[iv]*Red[i1][i2][iv];
   
        endfor
       
        // This the flux or O<-> as a function of the electrode reaction.
        JO[i1+1][iv] = (kf[i1+1][iv]*Ox[i1+1][2][iv] - kr[i1+1]*Red[i1+1][2][iv]) / (1 + Dx[iv]/D*(kf[i1+1][iv] + kr[i1+1][iv]));
        // This changes the concentration of O and R based on the flux.
        Ox[i1+1][1][iv] = Ox[i1+1][2][iv] - JO[i1+1][iv]*(Dx[iv]/D);
        // The km term in the red equation is in case R is degrading via an EC mechanism.
    Red[i1+1][1][iv] = Red[i1+1][2][iv] + JO[i1+1][iv]*(Dx[iv]/D) - km[iv]*Red[i1+1][1][iv];
        // Calculates current based on the flux - converted to mA by multiplying by 1000.
        Current[i1] = -n*F*A*C*JO[i1+1][iv];
       
    endfor
   
    // This is a quick-fix for naming each the current at each scan rate as a function of scan rate to make 1-d waves.
    String current_name = "current"+num2str(iv);   
    // This copies the "current" waves to its permanant wave "current[scan rate number]"
    Duplicate/O Current, $current_name
   
endfor

// Deleting all the extra waves created during the simulation.
killwaves Current, Dt, Dx, Enorm, E1, E2, JO, k, kf, km, kr, ktk, Ox, Red, tk, t_;

End