[lout,kprint,x,y,f,fx,fy,xe,ye,fe,de,fe2,fail]=evpcck(lout,kprint,x,y,f,fx,fy,xe,ye,fe,de,fe2,fail);

function [lout,kprint,x,y,f,fx,fy,xe,ye,fe,de,fe2,fail]=evpcck(lout,kprint,x,y,f,fx,fy,xe,ye,fe,de,fe2,fail);
persistent ax ay dermax derr dtrue dx faild faile failoc fdiff fdifmx fermax ferr firstCall ftrue i ier2 ierr inc j k machep ne nerr nmax nx ny pdermx pdifmx pfermx skip tol zero ; if isempty(firstCall),firstCall=1;end; 

f_orig=f;f_shape=[10,10];f=reshape([f_orig(1:min(prod(f_shape),numel(f_orig))),zeros(1,max(0,prod(f_shape)-numel(f_orig)))],f_shape);
fx_orig=fx;fx_shape=[10,10];fx=reshape([fx_orig(1:min(prod(fx_shape),numel(fx_orig))),zeros(1,max(0,prod(fx_shape)-numel(fx_orig)))],fx_shape);
fy_orig=fy;fy_shape=[10,10];fy=reshape([fy_orig(1:min(prod(fy_shape),numel(fy_orig))),zeros(1,max(0,prod(fy_shape)-numel(fy_orig)))],fy_shape);
if isempty(i), i=0; end;
if isempty(ier2), ier2=0; end;
if isempty(ierr), ierr=0; end;
if isempty(inc), inc=0; end;
if isempty(j), j=0; end;
if isempty(k), k=0; end;
if isempty(ne), ne=0; end;
if isempty(nerr), nerr=0; end;
if isempty(nmax), nmax=0; end;
if isempty(nx), nx=0; end;
if isempty(ny), ny=0; end;
if isempty(faild), faild=false; end;
if isempty(faile), faile=false; end;
if isempty(failoc), failoc=false; end;
if isempty(skip), skip=false; end;
if isempty(dermax), dermax=0; end;
if isempty(derr), derr=0; end;
if isempty(dtrue), dtrue=0; end;
if isempty(dx), dx=0; end;
if isempty(fdiff), fdiff=0; end;
if isempty(fdifmx), fdifmx=0; end;
if isempty(fermax), fermax=0; end;
if isempty(ferr), ferr=0; end;
if isempty(ftrue), ftrue=0; end;
if isempty(machep), machep=0; end;
if isempty(tol), tol=0; end;
if isempty(pdermx), pdermx=0; end;
if isempty(pdifmx), pdifmx=0; end;
if isempty(pfermx), pfermx=0; end;
if isempty(zero), zero=0; end;
if isempty(ax), ax=0; end;
if isempty(ay), ay=0; end;
% fcn= @(ax,ay)  ax*(ay*ay)*(ax*ax+1.0e0);real :: fcn ;
% dfdx= @(ax,ay) (ay*ay)*(3.0e0*ax*ax+1.0e0);real :: dfdx ;
% dfdy= @(ax,ay)  2.0e0*ax*ay*(ax*ax+1.0e0);real :: dfdy;
if firstCall,   nmax=[10];  end;
if firstCall,  nx=[4];  end;
if firstCall,  ny=[6];  end;
if firstCall,   ne=[51];  end;
if firstCall,   zero=[0.0e0];  end;
firstCall=0;
fcn= @(ax,ay)  ax.*(ay.*ay).*(ax.*ax+1.0e0);
dfdx= @(ax,ay) (ay.*ay).*(3.0e0.*ax.*ax+1.0e0);
dfdy= @(ax,ay)  2.0e0.*ax.*ay.*(ax.*ax+1.0e0);
[machep ]=r1mach(4);
tol = 10.0e0.*machep;
fail = false;
for i = 1 : nx - 1;
x(i) = 0.25e0.*i;
end; i = fix(nx - 1+1);
x(nx) = 1.0e0;
for j = 1 : ny;
y(j) = 0.5e0.*j - 1.25e0;
for i = 1 : nx;
fcn= @(ax,ay)  ax.*(ay.*ay).*(ax.*ax+1.0e0);
f(i,j) = fcn(x(i),y(j));
dfdx= @(ax,ay) (ay.*ay).*(3.0e0.*ax.*ax+1.0e0);
fx(i,j) = dfdx(x(i),y(j));
dfdy= @(ax,ay)  2.0e0.*ax.*ay.*(ax.*ax+1.0e0);
fy(i,j) = dfdy(x(i),y(j));
end; i = fix(nx+1);
end; j = fix(ny+1);
dx = 1.0e0./(ne-1);
for k = 1 : ne - 1;
xe(k) = dx.*(k-1);
ye(k) = 4.0e0.*xe(k) - 2.0e0;
end; k = fix(ne - 1+1);
xe(ne) = 1.0e0;
ye(ne) = 2.0e0;
if( kprint>=3 )
writef(lout,['1', '\n ' , '\n ' ,repmat(' ',1,10),'TEST PCHFE AND PCHFD' ' \n']);
end;
%format ('1'//10X,'TEST PCHFE AND PCHFD');
if( kprint>=2 )
writef(lout,[ '\n ' , '\n ' ,repmat(' ',1,10),'EVPCCK RESULTS', '\n ' ,repmat(' ',1,10),'--------------' ' \n']);
end;
%format [/10X,'EVPCCK RESULTS'/10X,'--------------');
nerr = 0;
inc = 1;
skip = false;
for j = 1 : ny;
[nx,x,f(sub2ind(size(f),1,j):end),fx(sub2ind(size(fx),1,j):end),inc,skip,ne,xe,fe,de,ierr]=pchfd(nx,x,f(sub2ind(size(f),1,j):end),fx(sub2ind(size(fx),1,j):end),inc,skip,ne,xe,fe,de,ierr);
if( kprint>=3 )
writef(lout,[ '\n ' , '\n ' ,repmat(' ',1,20),'PCHFD INCREMENT TEST -- INCFD = ','%2i', '\n ' ,repmat(' ',1,15),'ON ','%1s','-LINE ','%2i',',  ','%1s',' =','%8.4f','  --  IERR =','%3i' ' \n'], inc , 'J' , j , 'Y' , y(j) ,ierr);
end;
if( ierr<0 )
failoc = true;
if( kprint>=2 )
writef(lout,[ '\n ' , '\n ' ,'  ERROR  PCHFD RETURNED IERR =','%5i', '\n ' , '\n '  ' \n'], ierr);
end;
else;
if( kprint>3 )
writef(lout,[ '\n ' ,repmat(' ',1,3),'%1s','E',repmat(' ',1,10),'F',repmat(' ',1,8),'FE',repmat(' ',1,9),'DIFF',repmat(' ',1,13),'D',repmat(' ',1,8),'DE',repmat(' ',1,9),'DIFF' ' \n'], 'X');
end;
[nx,x,f(sub2ind(size(f),1,j):end),fx(sub2ind(size(fx),1,j):end),inc,skip,ne,xe,fe2,ier2]=pchfe(nx,x,f(sub2ind(size(f),1,j):end),fx(sub2ind(size(fx),1,j):end),inc,skip,ne,xe,fe2,ier2);
for k = 1 : ne;
fcn= @(ax,ay)  ax.*(ay.*ay).*(ax.*ax+1.0e0);
ftrue = fcn(xe(k),y(j));
ferr = fe(k) - ftrue;
dfdx= @(ax,ay) (ay.*ay).*(3.0e0.*ax.*ax+1.0e0);
dtrue = dfdx(xe(k),y(j));
derr = de(k) - dtrue;
if( kprint>3 )
writef(lout,['%7.2f',repmat([repmat(' ',1,2),repmat('%10.5f',1,2),repmat('%f',1,1),'%15.5f',repmat('%f',1,0)] ,1,2) ' \n'], xe(k) , ftrue , fe(k) ,ferr , dtrue , de(k) , derr);
end;
if( k==1 )
fermax = abs(ferr);
pfermx = xe(1);
dermax = abs(derr);
pdermx = xe(1);
fdifmx = abs(fe2(1)-fe(1));
pdifmx = xe(1);
else;
ferr = abs(ferr);
if( ferr>fermax )
fermax = ferr;
pfermx = xe(k);
end;
derr = abs(derr);
if( derr>dermax )
dermax = derr;
pdermx = xe(k);
end;
fdiff = abs(fe2(k)-fe(k));
if( fdiff>fdifmx )
fdifmx = fdiff;
pdifmx = xe(k);
end;
end;
end; k = fix(ne+1);
faild =(fermax>tol) |(dermax>tol);
faile = fdifmx~=zero;
failoc = faild | faile |(ierr~=13) |(ier2~=ierr);
if( failoc && (kprint>=2) )
writef(lout,[ '\n ' ,'  PCHFD AND/OR PCHFE FAILED ON ','%1s','-LINE ','%1i',',  ','%1s',' =','%8.4f' ' \n'], 'J' , j ,'Y' , y(j));
end;
if((kprint>=3) ||(faild &&(kprint==2)) )
writef(lout,[ '\n ' ,repmat(' ',1,17),'  MAXIMUM ERROR IN FUNCTION =',repmat('%f',1,1),repmat('%f',1,1),'%13.5f',repmat('%f',1,0),' (AT','%6.2f','),', '\n ' ,repmat(' ',1,31),'IN DERIVATIVE =',repmat('%f',1,1),'%13.5f',repmat('%f',1,0),' (AT','%6.2f',').' ' \n'], fermax , pfermx , dermax , pdermx);
end;
if( faild &&(kprint>=2) )
writef(lout,['  *** BOTH SHOULD BE  <=  TOL =',repmat('%f',1,1),'%12.5f',' ***' ' \n'], tol);
end;
if((kprint>=3) ||(faile &&(kprint==2)) )
writef(lout,['  MAXIMUM DIFFERENCE BETWEEN PCHFE AND PCHFD =',repmat('%f',1,1),'%13.5f',repmat('%f',1,0),' (AT','%6.2f',').' ' \n'], fdifmx , pdifmx);
end;
if( (ierr~=13) && (kprint>=2) )
writef(lout,[ '\n ' ,'  PCHF','%1s',' RETURNED IERR = ','%2i',' INSTEAD OF ','%2i' ' \n'], 'D' ,ierr , 13);
end;
if( (ier2~=ierr) && (kprint>=2) )
writef(lout,[ '\n ' ,'  PCHF','%1s',' RETURNED IERR = ','%2i',' INSTEAD OF ','%2i' ' \n'], 'E' ,ier2 , ierr);
end;
end;
if( failoc )
nerr = fix(nerr + 1);
end;
fail = fail | failoc;
end; j = fix(ny+1);
if( kprint>=2 )
if( nerr>0 )
writef(lout,[ '\n ' , '\n ' ,'  ERROR  PCHFD AND/OR PCHFE FAILED ON','%2i',repmat(' ',1,1),'%1s','-LINES.', '\n ' , '\n '  ' \n'], nerr , 'J');
else;
writef(lout,[ '\n ' ,' PCHFD AND PCHFE OK ON ','%1s','-LINES.' ' \n'], 'J');
end;
end;
nerr = 0;
inc = fix(nmax);
skip = false;
for i = 1 : nx;
[ny,y,f(sub2ind(size(f),i,1):end),fy(sub2ind(size(fy),i,1):end),inc,skip,ne,ye,fe,de,ierr]=pchfd(ny,y,f(sub2ind(size(f),i,1):end),fy(sub2ind(size(fy),i,1):end),inc,skip,ne,ye,fe,de,ierr);
if( kprint>=3 )
writef(lout,[ '\n ' , '\n ' ,repmat(' ',1,20),'PCHFD INCREMENT TEST -- INCFD = ','%2i', '\n ' ,repmat(' ',1,15),'ON ','%1s','-LINE ','%2i',',  ','%1s',' =','%8.4f','  --  IERR =','%3i' ' \n'], inc , 'I' , i , 'X' , x(i) ,ierr);
end;
if( ierr<0 )
failoc = true;
if( kprint>=2 )
writef(lout,[ '\n ' , '\n ' ,'  ERROR  PCHFD RETURNED IERR =','%5i', '\n ' , '\n '  ' \n'], ierr);
end;
else;
if( kprint>3 )
writef(lout,[ '\n ' ,repmat(' ',1,3),'%1s','E',repmat(' ',1,10),'F',repmat(' ',1,8),'FE',repmat(' ',1,9),'DIFF',repmat(' ',1,13),'D',repmat(' ',1,8),'DE',repmat(' ',1,9),'DIFF' ' \n'], 'Y');
end;
[ny,y,f(sub2ind(size(f),i,1):end),fy(sub2ind(size(fy),i,1):end),inc,skip,ne,ye,fe2,ier2]=pchfe(ny,y,f(sub2ind(size(f),i,1):end),fy(sub2ind(size(fy),i,1):end),inc,skip,ne,ye,fe2,ier2);
for k = 1 : ne;
fcn= @(ax,ay)  ax.*(ay.*ay).*(ax.*ax+1.0e0);
ftrue = fcn(x(i),ye(k));
ferr = fe(k) - ftrue;
dfdy= @(ax,ay)  2.0e0.*ax.*ay.*(ax.*ax+1.0e0);
dtrue = dfdy(x(i),ye(k));
derr = de(k) - dtrue;
if( kprint>3 )
writef(lout,['%7.2f',repmat([repmat(' ',1,2),repmat('%10.5f',1,2),repmat('%f',1,1),'%15.5f',repmat('%f',1,0)] ,1,2) ' \n'], ye(k) , ftrue , fe(k) ,ferr , dtrue , de(k) , derr);
end;
if( k==1 )
fermax = abs(ferr);
pfermx = ye(1);
dermax = abs(derr);
pdermx = ye(1);
fdifmx = abs(fe2(1)-fe(1));
pdifmx = ye(1);
else;
ferr = abs(ferr);
if( ferr>fermax )
fermax = ferr;
pfermx = ye(k);
end;
derr = abs(derr);
if( derr>dermax )
dermax = derr;
pdermx = ye(k);
end;
fdiff = abs(fe2(k)-fe(k));
if( fdiff>fdifmx )
fdifmx = fdiff;
pdifmx = ye(k);
end;
end;
end; k = fix(ne+1);
faild =(fermax>tol) |(dermax>tol);
faile = fdifmx~=zero;
failoc = faild | faile |(ierr~=20) |(ier2~=ierr);
if( failoc && (kprint>=2) )
writef(lout,[ '\n ' ,'  PCHFD AND/OR PCHFE FAILED ON ','%1s','-LINE ','%1i',',  ','%1s',' =','%8.4f' ' \n'], 'I' , i ,'X' , x(i));
end;
if((kprint>=3) ||(faild &&(kprint==2)) )
writef(lout,[ '\n ' ,repmat(' ',1,17),'  MAXIMUM ERROR IN FUNCTION =',repmat('%f',1,1),repmat('%f',1,1),'%13.5f',repmat('%f',1,0),' (AT','%6.2f','),', '\n ' ,repmat(' ',1,31),'IN DERIVATIVE =',repmat('%f',1,1),'%13.5f',repmat('%f',1,0),' (AT','%6.2f',').' ' \n'], fermax , pfermx , dermax , pdermx);
end;
if( faild &&(kprint>=2) )
writef(lout,['  *** BOTH SHOULD BE  <=  TOL =',repmat('%f',1,1),'%12.5f',' ***' ' \n'], tol);
end;
if((kprint>=3) ||(faile &&(kprint==2)) )
writef(lout,['  MAXIMUM DIFFERENCE BETWEEN PCHFE AND PCHFD =',repmat('%f',1,1),'%13.5f',repmat('%f',1,0),' (AT','%6.2f',').' ' \n'], fdifmx , pdifmx);
end;
if( (ierr~=20) && (kprint>=2) )
writef(lout,[ '\n ' ,'  PCHF','%1s',' RETURNED IERR = ','%2i',' INSTEAD OF ','%2i' ' \n'], 'D' ,ierr , 20);
end;
if( (ier2~=ierr) && (kprint>=2) )
writef(lout,[ '\n ' ,'  PCHF','%1s',' RETURNED IERR = ','%2i',' INSTEAD OF ','%2i' ' \n'], 'E' ,ier2 , ierr);
end;
end;
if( failoc )
nerr = fix(nerr + 1);
end;
fail = fail | failoc;
end; i = fix(nx+1);
if( kprint>=2 )
if( nerr>0 )
writef(lout,[ '\n ' , '\n ' ,'  ERROR  PCHFD AND/OR PCHFE FAILED ON','%2i',repmat(' ',1,1),'%1s','-LINES.', '\n ' , '\n '  ' \n'], nerr , 'I');
else;
writef(lout,[ '\n ' ,' PCHFD AND PCHFE OK ON ','%1s','-LINES.' ' \n'], 'I');
end;
end;
f_orig(1:prod(f_shape))=f;f=f_orig;
fx_orig(1:prod(fx_shape))=fx;fx=fx_orig;
fy_orig(1:prod(fy_shape))=fy;fy=fy_orig;
return;
%format [/20X,'PCHFD INCREMENT TEST -- INCFD = ',i2/15X,'ON ',a1,'-LINE ',i2,',  ',a1,' =',f8.4,'  --  IERR =',i3);
%format [3X,a1,'E',10X,'F',8X,'FE',9X,'DIFF',13X,'D',8X,'DE',9X,'DIFF');
%format(f7.2,2(2x,2f10.5,1p,e15.5,0p));
%format ['  PCHFD AND/OR PCHFE FAILED ON ',a1,'-LINE ',i1,',  ',a1,' =',f8.4);
%format [17X,'  MAXIMUM ERROR IN FUNCTION =',1P,1P,e13.5,0P,' (AT',f6.2,'),'/31X,'IN DERIVATIVE =',1P,e13.5,0P,' (AT',f6.2,').');
%format ('  MAXIMUM DIFFERENCE BETWEEN PCHFE AND PCHFD =',1P,e13.5,0P,' (AT',f6.2,').');
%format ['  PCHF',a1,' RETURNED IERR = ',i2,' INSTEAD OF ',i2);
%format ('  *** BOTH SHOULD BE  <=  TOL =',1P,e12.5,' ***');
%format [/'  ERROR  PCHFD RETURNED IERR =',i5/];
%format [/'  ERROR  PCHFD AND/OR PCHFE FAILED ON',i2,1X,a1,'-LINES.'/];
%format [' PCHFD AND PCHFE OK ON ',a1,'-LINES.');
f_orig(1:prod(f_shape))=f;f=f_orig;
fx_orig(1:prod(fx_shape))=fx;fx=fx_orig;
fy_orig(1:prod(fy_shape))=fy;fy=fy_orig;
end %subroutine evpcck