diff --git a/+gemini3d/+efield/Efield_BCs.m b/+gemini3d/+efield/Efield_BCs.m
index 2efc4e8..5ef53f4 100644
--- a/+gemini3d/+efield/Efield_BCs.m
+++ b/+gemini3d/+efield/Efield_BCs.m
@@ -61,12 +61,12 @@ function Efield_BCs(p, xg)
   end
 end
 
-thetamin = min(xg.theta(:));
-thetamax = max(xg.theta(:));
-mlatmin = 90-thetamax*180/pi;
-mlatmax = 90-thetamin*180/pi;
-mlonmin = min(xg.phi(:))*180/pi;
-mlonmax = max(xg.phi(:))*180/pi;
+thetamin = min(xg.theta, 'all');
+thetamax = max(xg.theta, 'all');
+mlatmin = 90- rad2deg(thetamax);
+mlatmax = 90- rad2deg(thetamin);
+mlonmin = rad2deg(min(xg.phi, 'all'));
+mlonmax = rad2deg(max(xg.phi, 'all'));
 
 % add a 1% buff
 latbuf = 1/100 * (mlatmax-mlatmin);
diff --git a/+gemini3d/+grid/model2geocoords.m b/+gemini3d/+grid/model2geocoords.m
index 23b356a..025a751 100644
--- a/+gemini3d/+grid/model2geocoords.m
+++ b/+gemini3d/+grid/model2geocoords.m
@@ -5,7 +5,7 @@
   lalt (1,1) {mustBePositive,mustBeInteger} = 150
   llon (1,1) {mustBePositive,mustBeInteger} = 150
   llat (1,1) {mustBePositive,mustBeInteger} = 150
-  altlims (1,2) {mustBeReal} = [min(xg.alt, 'all')+0.0001, max(xg.alt, 'all')-0.0001]
+  altlims (1,2) {mustBeReal} = [min(xg.alt, [], 'all')+0.0001, max(xg.alt, [], 'all')-0.0001]
   glonlims (1,2) {mustBeReal} = [nan,nan]
   glatlims (1,2) {mustBeReal} = [nan,nan]
 end
@@ -36,10 +36,10 @@
 %end %if
 
 if any(isnan(glonlims))
-  glonlims=[min(double(glon(:))), max(double(glon(:)))];
+  glonlims=[min(glon, [], 'all'), max(glon, [], 'all')];
 end
 if any(isnan(glatlims))
-  glatlims=[min(double(glat(:))), max(double(glat(:)))];
+  glatlims=[min(glat, [], 'all'), max(glat, [], 'all')];
 end
 
 %% Define a regular mesh of a set number of points that encompasses the grid (or part of the grid)
@@ -50,8 +50,8 @@
 
 
 %% Identify the type of grid that we are using
-minh1=min(xg.h1(:));
-maxh1=max(xg.h1(:));
+minh1=min(xg.h1, [], 'all');
+maxh1=max(xg.h1, [], 'all');
 if (abs(minh1-1)>1e-4 || abs(maxh1-1)>1e-4)    %curvilinear grid
   flagcurv=1;
 else                                           %cartesian grid
diff --git a/+gemini3d/+grid/model2magUENcoords.m b/+gemini3d/+grid/model2magUENcoords.m
index ecc7a87..7b8bc22 100644
--- a/+gemini3d/+grid/model2magUENcoords.m
+++ b/+gemini3d/+grid/model2magUENcoords.m
@@ -23,9 +23,9 @@
     lz=150; lx=150; ly=150;
 end %if
 if (nargin<8)    %default to using grid limits if not given
-    zlims=[min(zUEN(:))+1,max(zUEN(:))-1];   %stay just inside given grid
-    xlims=[min(xUEN(:))+1,max(xUEN(:))-1];
-    ylims=[min(yUEN(:))+1,max(yUEN(:))-1];
+    zlims=[min(zUEN, 'all')+1,max(zUEN, 'all')-1];   %stay just inside given grid
+    xlims=[min(xUEN, 'all')+1,max(xUEN, 'all')-1];
+    ylims=[min(yUEN, 'all')+1,max(yUEN, 'all')-1];
 end %if
 
 
@@ -37,8 +37,8 @@
 
 
 %% Identify the type of grid that we are using
-minh1=min(xg.h1(:));
-maxh1=max(xg.h1(:));
+minh1=min(xg.h1, 'all');
+maxh1=max(xg.h1, 'all');
 if (abs(minh1-1)>1e-4 || abs(maxh1-1)>1e-4)    %curvilinear grid
     flagcurv=1;
 else                                           %cartesian grid
diff --git a/+gemini3d/+model/magcalc.m b/+gemini3d/+model/magcalc.m
index d8fd3c3..78b5b39 100644
--- a/+gemini3d/+model/magcalc.m
+++ b/+gemini3d/+model/magcalc.m
@@ -39,8 +39,8 @@ function magcalc(direc, dang, xg)
 
   %TABULATE THE SOURCE OR GRID CENTER LOCATION
   if isempty(cfg.sourcemlon)
-    thdist=mean(xg.theta(:));
-    phidist=mean(xg.phi(:));
+    thdist=mean(xg.theta, 'all');
+    phidist=mean(xg.phi, 'all');
   else
     thdist= pi/2 - deg2rad(cfg.sourcemlat);    %zenith angle of source location
     phidist= deg2rad(cfg.sourcemlon);
diff --git a/+gemini3d/+particles/grid.m b/+gemini3d/+particles/grid.m
index 21ee827..a66b59b 100644
--- a/+gemini3d/+particles/grid.m
+++ b/+gemini3d/+particles/grid.m
@@ -64,12 +64,12 @@
 pg.Qit = zeros(pg.llon, pg.llat, Nt);
 pg.E0it = nan(pg.llon, pg.llat, Nt);
 
-thetamin = min(xg.theta(:));
-thetamax = max(xg.theta(:));
-mlatmin = 90-thetamax*180/pi;
-mlatmax = 90-thetamin*180/pi;
-mlonmin = min(xg.phi(:))*180/pi;
-mlonmax = max(xg.phi(:))*180/pi;
+thetamin = min(xg.theta, 'all');
+thetamax = max(xg.theta, 'all');
+mlatmin = 90- rad2deg(thetamax);
+mlatmax = 90- rad2deg(thetamin);
+mlonmin = rad2deg(min(xg.phi, 'all'));
+mlonmax = rad2deg(max(xg.phi, 'all'));
 
 % add a 1% buffer
 latbuf = 1/100*(mlatmax-mlatmin);
diff --git a/+gemini3d/+plot/cart2d.m b/+gemini3d/+plot/cart2d.m
index a6902a1..e09ef8c 100644
--- a/+gemini3d/+plot/cart2d.m
+++ b/+gemini3d/+plot/cart2d.m
@@ -46,8 +46,8 @@ function cart2d(time,xg,parm,parmlbl,caxlims,sourceloc, h, cmap)
 
 
 %SIZE OF PLOT GRID THAT WE ARE INTERPOLATING ONTO
-meantheta=mean(xg.theta(:));
-%meanphi=mean(xg.phi(:));
+meantheta=mean(xg.theta, 'all');
+%meanphi=mean(xg.phi, 'all');
 y=-1*(xg.theta-meantheta);   %this is a mag colat. coordinate and is only used for defining grid in linspaces below, runs backward from north distance, hence the negative sign
 %x=(xg.phi-meanphi);       %mag. lon coordinate, pos. eastward
 x=xg.x2(inds2)/Re/sin(meantheta);
@@ -55,12 +55,12 @@ function cart2d(time,xg,parm,parmlbl,caxlims,sourceloc, h, cmap)
 lxp=500;
 lyp=500;
 lzp=500;
-minx=min(x(:));
-maxx=max(x(:));
-miny=min(y(:));
-maxy=max(y(:));
-minz=min(z(:));
-maxz=max(z(:));
+minx=min(x, 'all');
+maxx=max(x, 'all');
+miny=min(y, 'all');
+maxy=max(y, 'all');
+minz=min(z, 'all');
+maxz=max(z, 'all');
 xp=linspace(minx,maxx,lxp);     %eastward distance (rads.)
 yp=linspace(miny,maxy,lyp);     %should be interpreted as northward distance (in rads.).  Irrespective of ordering of xg.theta, this will be monotonic increasing!!!
 zp=linspace(minz,maxz,lzp)';     %altitude (meters)
diff --git a/+gemini3d/+plot/cart3d_long.m b/+gemini3d/+plot/cart3d_long.m
index 3ef6c31..ec2f078 100644
--- a/+gemini3d/+plot/cart3d_long.m
+++ b/+gemini3d/+plot/cart3d_long.m
@@ -46,20 +46,20 @@ function cart3d_long(time, xg, parm, parmlbl, caxlims, sourceloc, hf, cmap)
 plotparams.altref=110;
 
 %% SIZE OF PLOT GRID THAT WE ARE INTERPOLATING ONTO
-meantheta=mean(xg.theta(:));
-meanphi=mean(xg.phi(:));
+meantheta=mean(xg.theta, 'all');
+meanphi=mean(xg.phi, 'all');
 y=-1*(xg.theta-meantheta);   %this is a mag colat. coordinate and is only used for defining grid in linspaces below, runs backward from north distance, hence the negative sign
 x=(xg.phi-meanphi);       %mag. lon coordinate, pos. eastward
 z=xg.alt/1e3;
 lxp=500;
 lyp=500;
 lzp=500;
-minx=min(x(:));
-maxx=max(x(:));
-miny=min(y(:));
-maxy=max(y(:));
-minz=min(z(:));
-maxz=max(z(:));
+minx=min(x, 'all');
+maxx=max(x, 'all');
+miny=min(y, 'all');
+maxy=max(y, 'all');
+minz=min(z, 'all');
+maxz=max(z, 'all');
 xp=linspace(minx,maxx,lxp);     %eastward distance (rads.)
 yp=linspace(miny,maxy,lyp);     %should be interpreted as northward distance (in rads.).  Irrespective of ordering of xg.theta, this will be monotonic increasing!!!
 zp=linspace(minz,maxz,lzp)';     %altitude (meters)
@@ -68,23 +68,23 @@ function cart3d_long(time, xg, parm, parmlbl, caxlims, sourceloc, hf, cmap)
 %ix1s=floor(lx1/2):lx1;    %only valide for a grid which is symmetric aboutu magnetic equator... (I think)
 ix1s=find(xg.x1(inds1)>=0);    %works for asymmetric grids
 minz=0;
-maxz=max(xg.alt(:));
+maxz=max(xg.alt, 'all');
 [tmp,ix1]=min(abs(xg.alt(ix1s,1,1)-maxz*1e3));
 ix1=ix1s(ix1);
 thetavals=xg.theta(ix1:lx1,:,:);
-meantheta=mean(thetavals(:));
+meantheta=mean(thetavals, 'all');
 phivals=xg.phi(ix1:lx1,:,:);
-meanphi=mean(phivals(:));
+meanphi=mean(phivals, 'all');
 x=(thetavals-meantheta);      %this is a mag colat. coordinate and is only used for defining grid in linspaces below and the parametric surfaces in the plots
 y=(phivals-meanphi);          %mag. lon coordinate
 z=xg.alt(ix1:lx1,:,:)/1e3;    %altitude
 lxp=500;
 lyp=500;
 lzp=500;
-minx=min(x(:));
-maxx=max(x(:));%+0.5*(max(x(:))-min(x(:)));
-miny=min(y(:));
-maxy=max(y(:));
+minx=min(x, 'all');
+maxx=max(x, 'all');%+0.5*(max(x, 'all')-min(x, 'all'));
+miny=min(y, 'all');
+maxy=max(y, 'all');
 xp=linspace(minx,maxx,lxp);
 yp=linspace(miny,maxy,lyp);
 zp=linspace(minz,maxz,lzp)';
diff --git a/+gemini3d/+plot/cart3d_long_enu.m b/+gemini3d/+plot/cart3d_long_enu.m
index 1c4590e..7ee1bdc 100644
--- a/+gemini3d/+plot/cart3d_long_enu.m
+++ b/+gemini3d/+plot/cart3d_long_enu.m
@@ -46,8 +46,8 @@ function cart3d_long_enu(time,xg,parm,parmlbl,caxlims,sourceloc,hf,cmap)
 plotparams.altref=300;
 
 %% SIZE OF PLOT GRID THAT WE ARE INTERPOLATING ONTO
-meantheta=mean(xg.theta(:));
-%meanphi=mean(xg.phi(:));
+meantheta=mean(xg.theta, 'all');
+%meanphi=mean(xg.phi, 'all');
 y=-1*(xg.theta-meantheta);   %this is a mag colat. coordinate and is only used for defining grid in linspaces below, runs backward from north distance, hence the negative sign
 %x=(xg.phi-meanphi);       %mag. lon coordinate, pos. eastward
 x=xg.x2(inds2)/Re/sin(meantheta);
@@ -55,12 +55,12 @@ function cart3d_long_enu(time,xg,parm,parmlbl,caxlims,sourceloc,hf,cmap)
 lxp=500;
 lyp=500;
 lzp=500;
-minx=min(x(:));
-maxx=max(x(:));
-miny=min(y(:));
-maxy=max(y(:));
-minz=min(z(:));
-maxz=max(z(:));
+minx=min(x, 'all');
+maxx=max(x, 'all');
+miny=min(y, 'all');
+maxy=max(y, 'all');
+minz=min(z, 'all');
+maxz=max(z, 'all');
 xp=linspace(minx,maxx,lxp);     %eastward distance (rads.)
 yp=linspace(miny,maxy,lyp);     %should be interpreted as northward distance (in rads.).  Irrespective of ordering of xg.theta, this will be monotonic increasing!!!
 zp=linspace(minz,maxz,lzp)';     %altitude (kilometers)
@@ -69,23 +69,23 @@ function cart3d_long_enu(time,xg,parm,parmlbl,caxlims,sourceloc,hf,cmap)
 %ix1s=floor(lx1/2):lx1;    %only valide for a grid which is symmetric aboutu magnetic equator... (I think)
 ix1s=find(xg.x1(inds1)>=0);    %works for asymmetric grids
 minz=0;
-maxz=max(xg.alt(:));
+maxz=max(xg.alt, 'all');
 [tmp,ix1]=min(abs(xg.alt(ix1s,1,1)-maxz*1e3));
 ix1=ix1s(ix1);
 thetavals=xg.theta(ix1:lx1,:,:);
-meantheta=mean(thetavals(:));
+meantheta=mean(thetavals, 'all');
 phivals=xg.phi(ix1:lx1,:,:);
-meanphi=mean(phivals(:));
+meanphi=mean(phivals, 'all');
 x=(thetavals-meantheta);      %this is a mag colat. coordinate and is only used for defining grid in linspaces below and the parametric surfaces in the plots
 y=(phivals-meanphi);          %mag. lon coordinate
 z=xg.alt(ix1:lx1,:,:)/1e3;    %altitude
 lxp=500;
 lyp=500;
 lzp=500;
-minx=min(x(:));
-maxx=max(x(:));%+0.5*(max(x(:))-min(x(:)));
-miny=min(y(:));
-maxy=max(y(:));
+minx=min(x, 'all');
+maxx=max(x, 'all');%+0.5*(max(x, 'all')-min(x, 'all'));
+miny=min(y, 'all');
+maxy=max(y, 'all');
 xp=linspace(minx,maxx,lxp);
 yp=linspace(miny,maxy,lyp);
 zp=linspace(minz,maxz,lzp)';
diff --git a/+gemini3d/+plot/curv2d.m b/+gemini3d/+plot/curv2d.m
index 4ce2489..d0c58b4 100644
--- a/+gemini3d/+plot/curv2d.m
+++ b/+gemini3d/+plot/curv2d.m
@@ -26,7 +26,7 @@ function curv2d(time,xg,parm,parmlbl,caxlims, sourceloc, h, cmap)
   sourcemlat=sourceloc(1);
   %sourcemlon=sourceloc(2);
 else
-  sourcemlat=mean(90-xg.theta(:)*180/pi);   %just take avg. over the grid
+  sourcemlat=mean(90-rad2deg(xg.theta), 'all');   %just take avg. over the grid
 end %if
 
 
@@ -45,8 +45,8 @@ function curv2d(time,xg,parm,parmlbl,caxlims, sourceloc, h, cmap)
 
 
 %SIZE OF PLOT GRID THAT WE ARE INTERPOLATING ONTO
-meantheta=mean(xg.theta(:));
-%meanphi=mean(xg.phi(:));
+meantheta=mean(xg.theta, 'all');
+%meanphi=mean(xg.phi, 'all');
 %meanphi=xg.x3(inds3(1));     %dont' forget that x3 has ghost cells!!!
 x=(xg.theta-meantheta);   %this is a mag colat. coordinate and is only used for defining grid in linspaces below
 %y=(xg.phi-meanphi);       %mag. lon coordinate
@@ -54,13 +54,13 @@ function curv2d(time,xg,parm,parmlbl,caxlims, sourceloc, h, cmap)
 lxp=1500;
 %lyp=500;
 lzp=1500;
-minx=min(x(:));
-maxx=max(x(:));
-%miny=min(y(:));
-%maxy=max(y(:));
-%minz=min(z(:));
+minx=min(x, 'all');
+maxx=max(x, 'all');
+%miny=min(y, 'all');
+%maxy=max(y, 'all');
+%minz=min(z, 'all');
 minz=0;
-maxz=max(z(:));
+maxz=max(z, 'all');
 xp=linspace(minx,maxx,lxp);
 %yp=linspace(miny,maxy,lyp);
 zp=linspace(minz,maxz,lzp)';
@@ -69,23 +69,23 @@ function curv2d(time,xg,parm,parmlbl,caxlims, sourceloc, h, cmap)
 %ix1s=floor(lx1/2):lx1;    %only valide for a grid which is symmetric aboutu magnetic equator... (I think)
 ix1s=find(xg.x1(inds1)>=0);    %works for asymmetric grids
 minz=0;
-maxz=max(xg.alt(:));
+maxz=max(xg.alt, 'all');
 [tmp,ix1]=min(abs(xg.alt(ix1s,1,1)-maxz*1e3));
 ix1=ix1s(ix1);
 thetavals=xg.theta(ix1:lx1,:,:);
-meantheta=mean(thetavals(:));
+meantheta=mean(thetavals, 'all');
 phivals=xg.phi(ix1:lx1,:,:);
-meanphi=mean(phivals(:));
+meanphi=mean(phivals, 'all');
 x=(thetavals-meantheta);      %this is a mag colat. coordinate and is only used for defining grid in linspaces below and the parametric surfaces in the plots
 y=(phivals-meanphi);          %mag. lon coordinate
 z=xg.alt(ix1:lx1,:,:)/1e3;    %altitude
 lxp=500;
 lyp=500;
 lzp=500;
-minx=min(x(:));
-maxx=max(x(:));%+0.5*(max(x(:))-min(x(:)));
-miny=min(y(:));
-maxy=max(y(:));
+minx=min(x, 'all');
+maxx=max(x, 'all');%+0.5*(max(x, 'all')-min(x, 'all'));
+miny=min(y, 'all');
+maxy=max(y, 'all');
 xp=linspace(minx,maxx,lxp);
 yp=linspace(miny,maxy,lyp);
 zp=linspace(minz,maxz,lzp)';
diff --git a/+gemini3d/+plot/curv3d_long.m b/+gemini3d/+plot/curv3d_long.m
index 33556ab..58f6923 100644
--- a/+gemini3d/+plot/curv3d_long.m
+++ b/+gemini3d/+plot/curv3d_long.m
@@ -66,8 +66,8 @@ function curv3d_long(time, xg, parm, parmlbl, caxlims, sourceloc, hf, cmap)
 plotparams.altref=375;
 
 %% SIZE OF PLOT GRID THAT WE ARE INTERPOLATING ONTO
-meantheta=mean(xg.theta(:));
-meanphi=mean(xg.phi(:));
+meantheta=mean(xg.theta, 'all');
+meanphi=mean(xg.phi, 'all');
 x=(xg.theta-meantheta);   %this is a mag colat. coordinate and is only used for defining grid in linspaces below
 y=(xg.phi-meanphi);       %mag. lon coordinate
 z=xg.alt/1e3;
@@ -75,13 +75,13 @@ function curv3d_long(time, xg, parm, parmlbl, caxlims, sourceloc, hf, cmap)
 lxp=1500;
 lyp=500;
 lzp=500;
-minx=min(x(:));
-maxx=max(x(:));
-miny=min(y(:));
-maxy=max(y(:));
-%minz=min(z(:));
+minx=min(x, 'all');
+maxx=max(x, 'all');
+miny=min(y, 'all');
+maxy=max(y, 'all');
+%minz=min(z, 'all');
 minz=-10e0;     %to give some space for the marker on th plots
-maxz=max(z(:));
+maxz=max(z, 'all');
 xp=linspace(minx,maxx,lxp);
 yp=linspace(miny,maxy,lyp);
 zp=linspace(minz,maxz,lzp)';
@@ -90,23 +90,23 @@ function curv3d_long(time, xg, parm, parmlbl, caxlims, sourceloc, hf, cmap)
 %ix1s=floor(lx1/2):lx1;    %only valide for a grid which is symmetric aboutu magnetic equator... (I think)
 ix1s=find(xg.x1(inds1)>=0);    %works for asymmetric grids
 minz=0;
-maxz=max(xg.alt(:));
+maxz=max(xg.alt, 'all');
 [tmp,ix1]=min(abs(xg.alt(ix1s,1,1)-maxz*1e3));
 ix1=ix1s(ix1);
 thetavals=xg.theta(ix1:lx1,:,:);
-meantheta=mean(thetavals(:));
+meantheta=mean(thetavals, 'all');
 phivals=xg.phi(ix1:lx1,:,:);
-meanphi=mean(phivals(:));
+meanphi=mean(phivals, 'all');
 x=(thetavals-meantheta);      %this is a mag colat. coordinate and is only used for defining grid in linspaces below and the parametric surfaces in the plots
 y=(phivals-meanphi);          %mag. lon coordinate
 z=xg.alt(ix1:lx1,:,:)/1e3;    %altitude
 lxp=500;
 lyp=500;
 lzp=500;
-minx=min(x(:));
-maxx=max(x(:));%+0.5*(max(x(:))-min(x(:)));
-miny=min(y(:));
-maxy=max(y(:));
+minx=min(x, 'all');
+maxx=max(x, 'all');%+0.5*(max(x, 'all')-min(x, 'all'));
+miny=min(y, 'all');
+maxy=max(y, 'all');
 xp=linspace(minx,maxx,lxp);
 yp=linspace(miny,maxy,lyp);
 zp=linspace(minz,maxz,lzp)';
@@ -161,7 +161,7 @@ function curv3d_long(time, xg, parm, parmlbl, caxlims, sourceloc, hf, cmap)
 if (flagsource)
   sourcetheta=pi/2-sourcemlat*pi/180;
 else
-  thetaref=mean(mean(xg.theta(1:floor(end/2),:,:)));
+  thetaref=mean(xg.theta(1:floor(end/2),:,:), 'all');
   sourcetheta=thetaref;
 end
 sourcex=sourcetheta-meantheta;
diff --git a/+gemini3d/+plot/diff.m b/+gemini3d/+plot/diff.m
index 79c6456..600bc93 100644
--- a/+gemini3d/+plot/diff.m
+++ b/+gemini3d/+plot/diff.m
@@ -80,7 +80,7 @@ function diff(A, B, name, time, newdir, refdir)
 plot(nexttile(t, 2), b)
 
 d = A - B;
-maxdiff = abs(max(d(:)));
+maxdiff = abs(max(d, 'all'));
 
 plot(nexttile(t, 3), d)
 
@@ -101,8 +101,8 @@ function diff(A, B, name, time, newdir, refdir)
   cmap = gemini3d.plot.bwr();
 end
 
-bmin = min(min(A(:)), min(B(:)));
-bmax = max(max(A(:)), max(B(:)));
+bmin = min(min(A, 'all'), min(B, 'all'));
+bmax = max(max(A, 'all'), max(B, 'all'));
 
 ax = nexttile(t, 1);
 hi = pcolor(ax, A);
@@ -132,9 +132,9 @@ function diff(A, B, name, time, newdir, refdir)
 %%
 ax = nexttile(t, 3);
 dAB = A - B;
-b = max(abs(min(dAB(:))), abs(max(dAB(:))));
+b = max(abs(min(dAB, 'all')), abs(max(dAB, 'all')));
 
-maxdiff = abs(max(dAB(:)));
+maxdiff = abs(max(dAB, 'all'));
 
 hi = pcolor(ax, dAB);
 set(hi, "EdgeColor", "none")
diff --git a/+gemini3d/+plot/frame.m b/+gemini3d/+plot/frame.m
index 1e368b3..ac555d3 100644
--- a/+gemini3d/+plot/frame.m
+++ b/+gemini3d/+plot/frame.m
@@ -69,41 +69,41 @@
 
 clm = opts.clim;
 if ~isfield(clm, "ne") && isfield(dat, 'ne')
-  clm.ne = [min(dat.ne(:)), max(dat.ne(:))];
+  clm.ne = [min(dat.ne, 'all'), max(dat.ne, 'all')];
 end
 if ~isfield(clm, "v1") && isfield(dat, 'v1')
-% v1mod=max(abs(v1(:)));
+% v1mod=max(abs(v1), 'all'));
   v1mod = 80;
   clm.v1 = [-v1mod, v1mod];
 end
 if ~isfield(clm, "Ti") && isfield(dat, 'Ti')
-  clm.Ti = [0, max(dat.Ti(:))];
+  clm.Ti = [0, max(dat.Ti, 'all')];
 end
 if ~isfield(clm, "Te") && isfield(dat, 'Te')
-  clm.Te = [0, max(dat.Te(:))];
+  clm.Te = [0, max(dat.Te, 'all')];
 end
 if ~isfield(clm, "J1") && isfield(dat, 'J1')
-  J1mod = max(abs(dat.J1(:)));
+  J1mod = max(abs(dat.J1), 'all');
   clm.J1 = [-J1mod, J1mod];
 end
 if ~isfield(clm, "v2") && isfield(dat, 'v2')
-  v2mod = max(abs(dat.v2(:)));
+  v2mod = max(abs(dat.v2), 'all');
   clm.v2 = [-v2mod, v2mod];
 end
 if ~isfield(clm, "v3") && isfield(dat, 'v3')
-  v3mod = max(abs(dat.v3(:)));
+  v3mod = max(abs(dat.v3), 'all');
   clm.v3 = [-v3mod, v3mod];
 end
 if ~isfield(clm, "J2") && isfield(dat, 'J2')
-  J2mod = max(abs(dat.J2(:)));
+  J2mod = max(abs(dat.J2), 'all');
   clm.J2 = [-J2mod, J2mod];
 end
 if ~isfield(clm, "J3") && isfield(dat, 'J3')
-  J3mod = max(abs(dat.J3(:)));
+  J3mod = max(abs(dat.J3), 'all');
   clm.J3 = [-J3mod, J3mod];
 end
 if ~isfield(clm, "Phitop") && isfield(dat, 'Phitop')
-  clm.Phitop = [min(dat.Phitop(:)), max(dat.Phitop(:))];
+  clm.Phitop = [min(dat.Phitop, 'all'), max(dat.Phitop, 'all')];
 end
 
 %% MAKE THE PLOTS
diff --git a/+gemini3d/+plot/grid2plotfun.m b/+gemini3d/+plot/grid2plotfun.m
index e9e3066..5d9f3e3 100644
--- a/+gemini3d/+plot/grid2plotfun.m
+++ b/+gemini3d/+plot/grid2plotfun.m
@@ -5,8 +5,8 @@
 
 %% DEFINE THE PLOTTING FUNCTION BASED ON THE TYPE OF GRID USED
 
-minh1=min(xg.h1(:));
-maxh1=max(xg.h1(:));
+minh1=min(xg.h1, 'all');
+maxh1=max(xg.h1, 'all');
 if (abs(minh1-1)>1e-4 || abs(maxh1-1)>1e-4)    %curvilinear grid
   if (xg.lx(2)>1 && xg.lx(3)>1)
     plotfun = @gemini3d.plot.curv3d_long;
diff --git a/+gemini3d/+plot/mag_map.m b/+gemini3d/+plot/mag_map.m
index 08ee4e2..8f87c41 100644
--- a/+gemini3d/+plot/mag_map.m
+++ b/+gemini3d/+plot/mag_map.m
@@ -72,8 +72,8 @@ function mag_map(direc)
 %INTERPOLATE TO HIGHER SPATIAL RESOLUTION FOR PLOTTING
 llonp=200;
 llatp=200;
-mlonp=linspace(min(mlon(:)),max(mlon(:)),llonp);
-mlatp=linspace(min(mlat(:)),max(mlat(:)),llatp);
+mlonp=linspace(min(mlon, 'all'),max(mlon, 'all'),llonp);
+mlatp=linspace(min(mlat, 'all'),max(mlat, 'all'),llatp);
 [MLONP,MLATP]=meshgrid(mlonp,mlatp);
 for it=1:length(cfg.times)
   param=interp2(cfg.mlon, cfg.mlat,squeeze(Brt(:,:,:,it)),MLONP,MLATP);
@@ -91,8 +91,8 @@ function mag_map(direc)
 
 %TABULATE THE SOURCE OR GRID CENTER LOCATION
 if isempty(cfg.sourcemlat)
-  thdist=mean(theta(:));
-  phidist=mean(phi(:));
+  thdist=mean(theta, 'all');
+  phidist=mean(phi, 'all');
   cfg.sourcemlat = 90 - rad2deg(thdist);
   cfg.sourcemlon = rad2deg(phidist);
 else
@@ -150,7 +150,7 @@ function plotBr(Brtp, cfg, ttxt)
 colormap(ax, gemini3d.plot.bwr());
 % set(ax,'FontSize',FS)
 tightmap
-caxlim=max(abs(param(:)));
+caxlim=max(abs(param), 'all');
 caxlim=max(caxlim,0.001);
 if verLessThan('matlab', '9.12')
   caxis(ax, [-caxlim,caxlim]) %#ok<CAXIS>
@@ -184,8 +184,8 @@ function plotBtheta(Bthetatp, cfg, ttxt)
 colormap(ax, gemini3d.plot.bwr())
 %set(ax,'FontSize',FS)
 tightmap
-caxlim=max(abs(param(:)));
-caxlim=max(caxlim,0.001);
+caxlim=max(abs(param), 'all');
+caxlim=max(caxlim, 0.001);
 if verLessThan('matlab', '9.12')
   caxis(ax, [-caxlim,caxlim]) %#ok<CAXIS>
 else
@@ -217,7 +217,7 @@ function plotBphi(Bphitp, cfg, ttxt)
 colormap(ax, gemini3d.plot.bwr())
 %set(ax,'FontSize',FS)
 tightmap
-caxlim=max(abs(param(:)));
+caxlim=max(abs(param), 'all');
 caxlim=max(caxlim,0.001);
 if verLessThan('matlab', '9.12')
   caxis(ax, [-caxlim,caxlim]) %#ok<CAXIS>
diff --git a/+gemini3d/+plot/mapgrid.m b/+gemini3d/+plot/mapgrid.m
index 8e9968a..53d484f 100644
--- a/+gemini3d/+plot/mapgrid.m
+++ b/+gemini3d/+plot/mapgrid.m
@@ -11,8 +11,8 @@
 % double for plot3m
 mlon = double(mlon);
 mlat = double(mlat);
-dmlon=max(mlon(:))-min(mlon(:));
-dmlat=max(mlat(:))-min(mlat(:));
+dmlon=max(mlon, 'all')-min(mlon, 'all');
+dmlat=max(mlat, 'all')-min(mlat, 'all');
 
 
 %% ORGANIZE INPUT STRUCTURE
@@ -41,8 +41,8 @@
     sourcemlonplot=sourcemlon;
   end
 else     %no "epicenter" to track just use mean grid locations
-  sourcemlon=mean(mlon(:));
-  sourcemlat=mean(mlat(:));
+  sourcemlon=mean(mlon, 'all');
+  sourcemlat=mean(mlat, 'all');
   sourcemlonplot=sourcemlon;
 end %if
 
@@ -70,9 +70,9 @@
 
 %% CONVERT INPUT GRID COORDINATES INTO MLAT,MLON,ALT
 % Re=6370e3;
-% dphi=max(xg.phi(:))-min(xg.phi(:));
+% dphi=max(xg.phi, 'all')-min(xg.phi, 'all');
 alt=xg.alt/1e3;
-altscale=max(alt(:));
+altscale=max(alt, 'all');
 alt=alt/altscale;
 % double for plot3m
 alt = double(alt);
@@ -264,7 +264,7 @@
     drho=rhomax-rhomin;
     xn=linspace(-1*drho,drho,lrho);
     rn=zn+6370;
-    dtheta=(max(xn(:))-min(xn(:)))/rn(1);
+    dtheta=(max(xn, 'all')-min(xn, 'all'))/rn(1);
     thetan=linspace(sourcetheta-dtheta/2,sourcetheta+dtheta/2,lrho);
     [THETAn,Rn]=meshgrid(thetan,rn);
 
@@ -291,7 +291,7 @@
 
       drho=rhomax-rhomin;                                                  %radius of circle, in kilometers, describing perp. directions of axisymmetric model
       xn=linspace(-1*drho,drho,lpts);                                       %N-S distance spanned by neutral model ("fake" number of grid points used here)
-      dthetan=(max(xn(:))-min(xn(:)))/rn(1);                               %equivalent theta coordinates of the neutral mesh (used in the plot of grid)
+      dthetan=(max(xn, 'all')-min(xn, 'all'))/rn(1);                               %equivalent theta coordinates of the neutral mesh (used in the plot of grid)
       thetan=linspace(sourcetheta-dthetan/2,sourcetheta+dthetan/2,lpts);    %theta coordinates of N-S distance specified
       phinhalf1=sourcephi+sqrt((dthetan/2)^2-(thetan-sourcetheta).^2);
       phinhalf2=sourcephi-sqrt((dthetan/2)^2-(thetan-sourcetheta).^2);
@@ -334,8 +334,8 @@
       yn=xn;              %arbitrarily let the y-extent be the same as x (our defs. of x,y left-handed here)
       rn=zn+6370;         %convert altitude to geocentric distance
 
-      dtheta=(max(xn(:))-min(xn(:)))/rn(1);    %equivalent theta coordinates of the neutral mesh
-      dphi=(max(yn(:))-min(yn(:)))/rn(1)/sin(sourcetheta);
+      dtheta=(max(xn, 'all')-min(xn, 'all'))/rn(1);    %equivalent theta coordinates of the neutral mesh
+      dphi=(max(yn, 'all')-min(yn, 'all'))/rn(1)/sin(sourcetheta);
       thetan=linspace(sourcetheta-dtheta/2,sourcetheta+dtheta/2,lpts);
       phin=linspace(sourcephi-dphi/2,sourcephi+dphi/2,lpts);
       [THETAn,PHIn,Rn]=meshgrid(thetan,phin,rn);
@@ -372,8 +372,8 @@
       yn=linspace(ymin,ymax,lpts);
       rn=zn+6370;         %convert altitude to geocentric distance
 
-      dtheta=(max(xn(:))-min(xn(:)))/rn(1);    %equivalent theta coordinates of the neutral mesh
-      dphi=(max(yn(:))-min(yn(:)))/rn(1)/sin(sourcetheta);
+      dtheta=(max(xn, 'all')-min(xn, 'all'))/rn(1);    %equivalent theta coordinates of the neutral mesh
+      dphi=(max(yn, 'all')-min(yn, 'all'))/rn(1)/sin(sourcetheta);
       thetan=linspace(sourcetheta-dtheta/2,sourcetheta+dtheta/2,lpts);
       phin=linspace(sourcephi-dphi/2,sourcephi+dphi/2,lpts);
       [THETAn,PHIn,Rn]=meshgrid(thetan,phin,rn);
diff --git a/example/TECcalc.m b/example/TECcalc.m
index 5a27c84..5636210 100644
--- a/example/TECcalc.m
+++ b/example/TECcalc.m
@@ -41,8 +41,8 @@
 
 %DEFINE A CENTER AND REGION OF INTEREST
 if (isempty(cfg.sourcemlat))    %in case this run didn't have a disturbance!
-  cfg.sourcemlat = (pi/2-mean(xg.theta(:)))*180/pi;
-  cfg.sourcemlon = mean(xg.phi(:))*180/pi;
+  cfg.sourcemlat = pi/2 - rad2deg(mean(xg.theta, 'all'));
+  cfg.sourcemlon = rad2deg(mean(xg.phi, 'all'));
 end
 thdist = pi/2 - deg2rad(cfg.sourcemlat);    %zenith angle of source location
 phidist = deg2rad(cfg.sourcemlon);
@@ -66,12 +66,12 @@
 rvals=xg.r(1:lh,:,:);
 thvals=xg.theta(1:lh,:,:);
 phivals=xg.phi(1:lh,:,:);
-rmin=min(rvals(:));
-rmax=max(rvals(:));
-thmin=min(thvals(:));
-thmax=max(thvals(:));
-phimin=min(phivals(:));
-phimax=max(phivals(:));
+rmin=min(rvals, [], 'all');
+rmax=max(rvals, [], 'all');
+thmin=min(thvals, [], 'all');
+thmax=max(thvals, [], 'all');
+phimin=min(phivals, [], 'all');
+phimax=max(phivals, [], 'all');
 
 theta=linspace(thmin,thmax,lth);
 r=linspace(rmin,rmax,lr)';
@@ -227,7 +227,7 @@
 %       set(gca,'FontSize',FS);
       axis xy;
       axis tight;
-      caxlim=max(max(abs(dvTEC(:,:,it))));
+      caxlim= max(abs(dvTEC(:,:,it)), [], 'all');
       caxlim=max(caxlim,0.01);
       if verLessThan('matlab', '9.12')
         caxis([-1*caxlim, caxlim]) %#ok<CAXIS>
@@ -264,9 +264,9 @@
   datetick;
   axis tight;
   if verLessThan('matlab', '9.12')
-    caxis([-max(max(abs(dvTEC(:,:)))), max(max(abs(dvTEC(:,:))))]) %#ok<CAXIS>
+    caxis([-max(abs(dvTEC), [], 'all'), max(abs(dvTEC), [], 'all')]) %#ok<CAXIS>
   else
-    clim([-max(max(abs(dvTEC(:,:)))), max(max(abs(dvTEC(:,:))))])
+    clim([-max(abs(dvTEC), [], 'all'), max(abs(dvTEC), [], 'all')])
   end
   c=colorbar;
 %   set(c,'FontSize',FS)