My main script contains following code:

%# Grid and model parameters
nModel=50;
nModel_want=1;
nI_grid1=5;
Nth=1;
nRow.Scale1=5;
nCol.Scale1=5;
nRow.Scale2=5^2;
nCol.Scale2=5^2;

theta = 90; % degrees
a_minor = 2; % range along minor direction
a_major = 5; % range along major direction
sill = var(reshape(Deff_matrix_NthModel,nCell.Scale1,1)); % variance of the coarse data matrix of size nRow.Scale1 X nCol.Scale1

%#  Covariance computation

% Scale 1
for ihRow = 1:nRow.Scale1
    for ihCol = 1:nCol.Scale1
        [cov.Scale1(ihRow,ihCol),heff.Scale1(ihRow,ihCol)] = general_CovModel(theta, ihCol, ihRow, a_minor, a_major, sill, 'Exp');
    end
end

% Scale 2
for ihRow = 1:nRow.Scale2
    for ihCol = 1:nCol.Scale2
        [cov.Scale2(ihRow,ihCol),heff.Scale2(ihRow,ihCol)] = general_CovModel(theta, ihCol/(nCol.Scale2/nCol.Scale1), ihRow/(nRow.Scale2/nRow.Scale1), a_minor, a_major, sill/(nRow.Scale2*nCol.Scale2), 'Exp');
    end
end


%# Scale-up of fine scale values by averaging
[covAvg.Scale2,var_covAvg.Scale2,varNorm_covAvg.Scale2] = general_AverageProperty(nRow.Scale2/nRow.Scale1,nCol.Scale2/nCol.Scale1,1,nRow.Scale1,nCol.Scale1,1,cov.Scale2,1);

I am using two functions, general_CovModel() and general_AverageProperty(), in my main script which are given as following:

function [cov,h_eff] = general_CovModel(theta, hx, hy, a_minor, a_major, sill, mod_type)
% mod_type should be in strings

angle_rad = theta*(pi/180); % theta in degrees, angle_rad in radians
R_theta = [sin(angle_rad) cos(angle_rad); -cos(angle_rad) sin(angle_rad)];
h = [hx; hy];
lambda = a_minor/a_major;
D_lambda = [lambda 0; 0 1];
h_2prime = D_lambda*R_theta*h;
h_eff = sqrt((h_2prime(1)^2)+(h_2prime(2)^2));

if strcmp(mod_type,'Sph')==1 || strcmp(mod_type,'sph') ==1
    if h_eff<=a
        cov = sill - sill.*(1.5*(h_eff/a_minor)-0.5*((h_eff/a_minor)^3));
    else
        cov = sill;
    end
elseif strcmp(mod_type,'Exp')==1 || strcmp(mod_type,'exp') ==1
    cov = sill-(sill.*(1-exp(-(3*h_eff)/a_minor)));
elseif strcmp(mod_type,'Gauss')==1 || strcmp(mod_type,'gauss') ==1
    cov = sill-(sill.*(1-exp(-((3*h_eff)^2/(a_minor^2)))));
end

and

function [PropertyAvg,variance_PropertyAvg,NormVariance_PropertyAvg]=...
    general_AverageProperty(blocksize_row,blocksize_col,blocksize_t,...
    nUpscaledRow,nUpscaledCol,nUpscaledT,PropertyArray,omega)
% This function computes average of a property and variance of that averaged
% property using power averaging

PropertyAvg=zeros(nUpscaledRow,nUpscaledCol,nUpscaledT);

%# Average of property
for k=1:nUpscaledT,
    for j=1:nUpscaledCol,
        for i=1:nUpscaledRow,
            sum=0;
            for a=1:blocksize_row,
                for b=1:blocksize_col,
                    for c=1:blocksize_t,
                        sum=sum+(PropertyArray((i-1)*blocksize_row+a,(j-1)*blocksize_col+b,(k-1)*blocksize_t+c).^omega); % add all the property values in 'blocksize_x','blocksize_y','blocksize_t' to one variable
                    end
                end
            end
            PropertyAvg(i,j,k)=(sum/(blocksize_row*blocksize_col*blocksize_t)).^(1/omega); % take average of the summed property
        end
    end
end

%# Variance of averageed property
variance_PropertyAvg=var(reshape(PropertyAvg,...
    nUpscaledRow*nUpscaledCol*nUpscaledT,1),1,1); 

%# Normalized variance of averageed property
NormVariance_PropertyAvg=variance_PropertyAvg./(var(reshape(...
    PropertyArray,numel(PropertyArray),1),1,1)); 

Question: Using Matlab, I would like to optimize covAvg.Scale2 such that it matches closely with cov.Scale1 by perturbing/varying any (or all) of the following variables

1) a_minor

2) a_major

3) theta

Thanks.