function [pca_eigenval]=PCAgenerator

disp('Multi-parameter variation in XPP from Matlab, with calculating principal components')
disp(' ')
disp('The parameters varied are wights w1, stimulation strength Kn, gsyn and iapp')
disp(' ')

% user specifications
odeFileName = 'tenStim3electrgap_LF.ode';

x1range = 0.3:0.1:0.3;
x2range = 0.0:2.0:60.0;
x3range = 0.5:0.1:1.4;
x4range = 4.0:1.0:8.0;

verboseTog = false; % two levels of notification verbosity to user

% initializations
x1_ix = 0;

lenx1 = length(x1range); 
lenx2 = length(x2range);
lenx3 = length(x3range);
lenx4 = length(x4range);

X2 = x2range;
X3 = x3range;
X4 = x4range;

numpoints = lenx1 * lenx2 * lenx3 * lenx4;

fprintf('Calculating %d points\n',numpoints)
pointnum = 0;
count = 0; % line print count, for backspacing purposes

newpars(1).name='w1';
newpars(1).type='PAR';
newpars(2).name='Kn';
newpars(2).type='PAR';
newpars(3).name='gsyn';
newpars(3).type='PAR';
newpars(4).name='iapp';
newpars(4).type='PAR';


for i=x1range
    x1_ix = x1_ix+1;
    pca_eigenval = cell(lenx2,lenx3,lenx4);
    x2_ix = 0;
    for j = x2range
        x2_ix = x2_ix+1;
        x3_ix = 0;
        for k = x3range
            x3_ix = x3_ix+1;
            x4_ix = 0;
            for l = x4range
                x4_ix = x4_ix+1;        
                pointnum = pointnum + 1;
                if ~verboseTog % then print current loop counter
                    bspstr = '';
                    if count > 0
                        for bix=1:count
                            bspstr = [bspstr '\b']; 
                        end
                    end
                    fprintf(bspstr);
                    count = fprintf('%d',pointnum);
                end
 newpars(1).val = i;
 newpars(2).val = j;
 newpars(3).val = k;
 newpars(4).val = l;
 if verboseTog
     fprintf(' **** Point %d: Variable pt = %.4f, c = %.4f\n',pointnum,newpars(1).val,newpars(2).val)
 end
 success = ChangeXPPodeFile(odeFileName,newpars); % type `help ChangeXPPodeFile` in Matlab
 if success==0
     disp('Change XPP ode file failed')
     return
 end
 success = RunXPP(odeFileName,'',true,'C:\xppall\xppaut'); % type `help RunXPP` in Matlab
 if success==0
     disp('Run XPP failed')
     return
 end
 format = '%f32%*f32%*f32%*f32%*f32%*f32%*f32%*f32%*f32%*f32%*f32%f32%f32%f32%*f32%*f32%*f32%f32%f32%f32%f32%f32%f32%f32%f32%f32%f32%*[^\n]';
 fid = fopen('output.dat', 'r'); % change this if your ODE file explicitly names a data file
 data = textscan(fid, format);
 fclose(fid);
  
 k_v5 = 0;
 k_v6 = 0;
 k_v7 = 0;
 k_v = 0;
 
 for m = 120001:139981 %count # of spikes during the last minute of time series
    if (data{2}(m-1) < -40) && (data{2}(m) >= -40)
        k_v5 = k_v5+1;
    end
    if (data{3}(m-1) < -40) && (data{3}(m) >= -40)
        k_v6 = k_v6+1;
    end
    if (data{4}(m-1) < -40) && (data{4}(m) >= -40)
        k_v7 = k_v7+1;
    end
    if (data{3}(m) < -200) || (data{3}(m) > 200)
        k_v=k_v+1;
    end
 end
 %compute PCA eigenvalues using the slow variable r
 if (min([k_v5 k_v6 k_v7]) > 10) && (max([k_v5 k_v6 k_v7]) < 100) && (k_v==0)
    pca_eigenval{x2_ix,x3_ix,x4_ix}=pca([data{5}(40001:139981) data{6}(40001:139981) data{7}(40001:139981) data{8}(40001:139981) data{9}(40001:139981) data{10}(40001:139981) data{11}(40001:139981) data{12}(40001:139981) data{13}(40001:139981) data{14}(40001:139981)]);
    %t1 = 6000:0.5:6990; left for plotting purposes
   % v5 = data{3}(12001:13981);
   % fOut1 = sprintf('Volt_tr_%s%1.1f_%s%2.1f_%s%1.1f_%s%1.1f.mat','w',i,'Kn',j,'gsyn',k,'iapp',l);
   % save(fOut1, 'v5')
   % clear v5
 end 
    clear data
             end
        end
    end
    fOut2 = sprintf('PCA_eigenval_%s%1.1f.mat','w',i);
    save(fOut2,'X4','X3','X2','pca_eigenval')
    clear pca_eigenval
end

fprintf('\n')
beep

return

end