diff --git a/tools/plot/matlab/build_subnet.m b/tools/plot/matlab/build_subnet.m
new file mode 100644
index 00000000..c176537f
--- /dev/null
+++ b/tools/plot/matlab/build_subnet.m
@@ -0,0 +1,51 @@
+function [zn_new,nname_new] = build_subnet(filename,min_flx,min_xmf,max_z)
+
+[zz, aa, xmf, time, ~, ~, timestep, ~, flx_end, flx] = read_ts_file(filename);
+cycle_limit = [1, size(time,2)];
+
+nn = aa - zz;
+zn = [zz, nn];
+
+% Maximum Z
+zn_zmax_new = zn(zz <= max_z,:);
+
+% Maximum of each species during integration
+xmf_max = max(xmf,[],2);
+zn_xmf_new = zn(xmf_max > min_xmf,:);
+
+% Max flux
+flx_max = zeros(size(flx,1),1);
+[~,iwrite] = max(abs(flx),[],2);
+for j = 1:size(flx,1)
+    flx_max(j) = flx(j,iwrite(j));
+end
+iflx_new = find(abs(flx_max) > min_flx*max(abs(flx_max)));
+znt = [zz(flx_end(:,1)), nn(flx_end(:,1))];
+znp = [zz(flx_end(:,2)), nn(flx_end(:,2))];
+znt_flx_new = znt(iflx_new,:);
+znp_flx_new = znp(iflx_new,:);
+
+% Construct new network: Z <= max_z and ( max(X) > min_xmf or flx_max > min_flx )
+zn_flx_new = unique( [znt_flx_new;znp_flx_new], 'rows' );
+zn_new = intersect(zn_zmax_new,union(zn_xmf_new,zn_flx_new,'rows'),'rows');
+
+element = build_element_symbol;
+
+% Build the species list (sunet)
+nname_new = cell(size(zn_new,1),1);
+for i = 1:size(zn_new,1)
+    label = sprintf('%s%d',lower(strtrim(element(zn_new(i,1)+1,:))),zn_new(i,1)+zn_new(i,2));
+    if strcmp(label,'n1')
+        label = 'n';
+    elseif strcmp(label,'h1')
+        label = 'p';
+    elseif strcmp(label,'h2')
+        label = 'd';
+    elseif strcmp(label,'h3')
+        label = 't';
+    end
+    nname_new(i) = {label};
+end
+
+end
+
diff --git a/tools/plot/matlab/compare_final.m b/tools/plot/matlab/compare_final.m
new file mode 100644
index 00000000..b6750ade
--- /dev/null
+++ b/tools/plot/matlab/compare_final.m
@@ -0,0 +1,45 @@
+function [aerr,rerr] = compare_final(ts1,ts2,varargin)
+%UNTITLED Summary of this function goes here
+%   Detailed explanation goes here
+
+if ( nargin > 2 )
+    verbose = varargin{1};
+else
+    verbose = false;
+end
+
+ny = ts1.ny;
+
+enuc1 = cumsum(ts1.edot{:} .* ts1.tdel{:});
+enuc2 = cumsum(ts2.edot{:} .* ts2.tdel{:});
+
+nuc_name = build_isotope_symbol(ts1.zz{:},ts1.aa{:});
+
+aerr = abs( ts2.xn{1}(:,end) - ts1.xn{1}(:,end) );
+rerr = aerr ./ ts1.xn{1}(:,end);
+
+[~,isort] = sort(rerr);
+aerr(ny+1) = norm(aerr,2);
+rerr(ny+1) = aerr(ny+1) ./ norm(ts1.xn{1}(:,end),2);
+
+aerr(ny+2) = abs(ts2.t9{1}(end) - ts1.t9{1}(end));
+rerr(ny+2) = aerr(ny+2) ./ ts1.t9{1}(end);
+
+aerr(ny+3) = abs(enuc2(end) - enuc1(end));
+rerr(ny+3) = aerr(ny+3) ./ enuc1(end);
+
+if ( verbose )
+    disp(sprintf("%5s %5s\t%10s\t%16s\t%16s\t%16s\t%16s",'i','isort','name','X1','X2','|dX|','|dX| / |X1|'));
+    fmt = "%5d %5d\t%10s\t%16.8e\t%16.8e\t%16.8e\t%16.8e";
+    for i = 1:ny
+        disp(sprintf(fmt, i, isort(i), nuc_name{isort(i)}, ts1.xn{1}(isort(i),end), ts2.xn{1}(isort(i),end), aerr(isort(i)), rerr(isort(i))));
+    end
+    fmt = "%5s %5s\t%10s\t%16s\t%16s\t%16.8e\t%16.8e";
+    disp(sprintf(fmt, '', '', '2-norm', '', '', aerr(ny+1), rerr(ny+1)));
+    fmt = "%5s %5s\t%10s\t%16.8e\t%16.8e\t%16.8e\t%16.8e";
+    disp(sprintf(fmt, '', '', 'T', ts1.t9{1}(end), ts2.t9{1}(end), aerr(ny+2), rerr(ny+2)));
+    disp(sprintf(fmt, '', '', 'E_nuc', enuc1(end), enuc2(end), aerr(ny+3), rerr(ny+3)));
+end
+
+end
+
diff --git a/tools/plot/matlab/draw_nz_network_grid.m b/tools/plot/matlab/draw_nz_network_grid.m
new file mode 100644
index 00000000..ec4063eb
--- /dev/null
+++ b/tools/plot/matlab/draw_nz_network_grid.m
@@ -0,0 +1,60 @@
+function [axis_h,group_h,line_h] = draw_nz_network_grid(data_directory)
+%--------------------------------------------------------------------------
+% [axis_id] = draw_nz_network(data_directory,pcolor,psize)
+% Scatter plot of the isotopes in a specified netwinv file.
+% Inputs>  data_directory: directory hosting netwinv file to be plotted.
+%          pcolor: point color
+%          psize: point size
+% Outputs: axis_id: handle of current axis
+%--------------------------------------------------------------------------
+% Open file
+  filename=strcat(data_directory,'/netwinv');
+  fileID=fopen(filename);
+  
+% Read network size
+  dataread=textscan(fileID,'%d',1);
+  ny=cell2mat(dataread);
+
+% Skip thermodata
+  dataread=textscan(fileID,'%d',1);
+  
+% Read nuclear names
+  dataread = textscan(fileID,'%s',ny);
+  nname=dataread{1};
+
+% Read nuclear data, skipping partition function  
+  file_form='%*s %f %f %f %f %f';
+  for i=1:ny        
+    dataread=textscan(fileID,file_form,1);
+    [aa(i),zz(i),nn(i),sp(i),be(i)]=dataread{1:5};
+    
+    dataread=textscan(fileID,'%f %f %f %f %f %f %f %f',3);
+
+  end
+  [minz,~] = min(zz);
+  [maxz,~] = max(zz);
+  line_h = [];
+  for iz=minz:maxz
+      [minn,~] = min(nn(zz==iz));
+      [maxn,~] = max(nn(zz==iz));
+      if length(minn) > 0 & length(maxn) > 0
+          line_h = [line_h; plot( [minn-0.5,maxn+0.5], [iz-0.5, iz-0.5], 'k' )];
+          line_h = [line_h; plot( [minn-0.5,maxn+0.5], [iz+0.5, iz+0.5], 'k' )];
+      end
+  end
+  [minn,~] = min(nn);
+  [maxn,~] = max(nn);
+  for in=minn:maxn
+      [minz,~] = min(zz(nn==in));
+      [maxz,~] = max(zz(nn==in));
+      if length(minz) > 0 & length(maxz) > 0
+          line_h = [line_h; plot( [in-0.5,in-0.5], [minz-0.5, maxz+0.5], 'k' )];
+          line_h = [line_h; plot( [in+0.5,in+0.5], [minz-0.5, maxz+0.5], 'k' )];
+      end
+  end
+  group_h = hggroup( 'Parent', gca );
+  set(get(get(group_h,'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
+  set(line_h,'Parent',group_h);
+  axis_h = gca;
+
+end
\ No newline at end of file
diff --git a/tools/plot/matlab/draw_nz_network_patch.m b/tools/plot/matlab/draw_nz_network_patch.m
new file mode 100644
index 00000000..bfbe3c35
--- /dev/null
+++ b/tools/plot/matlab/draw_nz_network_patch.m
@@ -0,0 +1,58 @@
+function [axis_h,patch_h] = draw_nz_network_patch(data_directory)
+%--------------------------------------------------------------------------
+% [axis_id] = draw_nz_network(data_directory,pcolor,psize)
+% Scatter plot of the isotopes in a specified netwinv file.
+% Inputs>  data_directory: directory hosting netwinv file to be plotted.
+%          pcolor: point color
+%          psize: point size
+% Outputs: axis_id: handle of current axis
+%--------------------------------------------------------------------------
+% Open file
+  filename=strcat(data_directory,'/netwinv');
+  fileID=fopen(filename);
+  
+% Read network size
+  dataread=textscan(fileID,'%d',1);
+  ny=cell2mat(dataread);
+
+% Skip thermodata
+  dataread=textscan(fileID,'%d',1);
+  
+% Read nuclear names
+  dataread = textscan(fileID,'%s',ny);
+  nname=dataread{1};
+
+% Read nuclear data, skipping partition function  
+  file_form='%*s %f %f %f %f %f';
+  for i=1:ny        
+    dataread=textscan(fileID,file_form,1);
+    [aa(i),zz(i),nn(i),sp(i),be(i)]=dataread{1:5};
+    
+    dataread=textscan(fileID,'%f %f %f %f %f %f %f %f',3);
+
+  end
+  [minz,~] = min(zz);
+  [maxz,~] = max(zz);
+  [minn,~] = min(nn(zz==minz));
+  [maxn,~] = max(nn(zz==minz));
+  xvertices = maxn+0.5:-1:minn-0.5;
+  yvertices = repmat(minz-0.5,1,maxn-minn+2);
+  for iz=minz:maxz
+      [minn,~] = min(nn(zz==iz));
+      xvertices = [xvertices,minn-0.5,minn-0.5];
+      yvertices = [yvertices,iz-0.5,iz+0.5];
+  end
+  [minn,~] = min(nn(zz==maxz));
+  [maxn,~] = max(nn(zz==maxz));
+  xvertices = [xvertices,minn-0.5:maxn+0.5];
+  yvertices = [yvertices,repmat(maxz+0.5,1,maxn-minn+2)];
+  for iz=maxz:-1:minz
+      [maxn,~] = max(nn(zz==iz));
+      xvertices = [xvertices,maxn+0.5,maxn+0.5];
+      yvertices = [yvertices,iz+0.5,iz-0.5];
+  end
+  patch_h = patch('XData',xvertices,'YData',yvertices,'FaceColor','none');
+  axis_h = gca;
+%   scatter(nn,zz,psize,'Marker','s','MarkerEdgeColor',pcolor,'MarkerFaceColor','none');
+
+end
\ No newline at end of file
diff --git a/tools/plot/matlab/read_nz_network.m b/tools/plot/matlab/read_nz_network.m
new file mode 100644
index 00000000..31cd1d9a
--- /dev/null
+++ b/tools/plot/matlab/read_nz_network.m
@@ -0,0 +1,40 @@
+function [nname,aa,zz,nn,sp,be] = read_nz_network(data_directory)
+%--------------------------------------------------------------------------
+% [axis_id] = draw_nz_network(data_directory,pcolor,psize)
+% Scatter plot of the isotopes in a specified netwinv file.
+% Inputs>  data_directory: directory hosting netwinv file to be plotted.
+%          pcolor: point color
+%          psize: point size
+% Outputs: axis_id: handle of current axis
+%--------------------------------------------------------------------------
+% Open file
+  filename=strcat(data_directory,'/netwinv');
+  fileID=fopen(filename);
+  
+% Read network size
+  dataread=textscan(fileID,'%d',1);
+  ny=cell2mat(dataread);
+
+% Skip thermodata
+  dataread=textscan(fileID,'%d',1);
+  
+% Read nuclear names
+  dataread = textscan(fileID,'%s',ny);
+  nname=dataread{1};
+
+% Read nuclear data, skipping partition function
+  aa = zeros(ny,1);
+  zz = zeros(ny,1);
+  nn = zeros(ny,1);
+  sp = zeros(ny,1);
+  be = zeros(ny,1);
+  file_form='%*s %f %f %f %f %f';
+  for i=1:ny        
+    dataread=textscan(fileID,file_form,1);
+    [aa(i),zz(i),nn(i),sp(i),be(i)]=dataread{1:5};
+    
+    dataread=textscan(fileID,'%f %f %f %f %f %f %f %f',3);
+
+  end
+
+end
\ No newline at end of file
diff --git a/tools/plot/matlab/read_sparse_ind.m b/tools/plot/matlab/read_sparse_ind.m
new file mode 100644
index 00000000..db17c3e4
--- /dev/null
+++ b/tools/plot/matlab/read_sparse_ind.m
@@ -0,0 +1,65 @@
+function [ny, nnz, ridx, cidx, pb, nrext, ns1, ns2, ns3, ns4 ] = read_sparse_ind( data_directory )
+%--------------------------------------------------------------------------
+%[zz, aa, xmf, time, temperature, density, timestep, edot, flx_end, flx] = read_ts_file( sparse_filename ) 
+% Reads XNet binary sparsity pattern file
+% Inputs>  sparse_filename: name of binary file 
+% Outputs< zz: proton number for each isotope in the network.
+%          aa: total nucleon number for each isotope in the network
+%          xmf: time series of mass fractions 
+%          time: discrete temporal evolution
+%          temperature:  time series of temperature
+%          density: time series of density
+%          timestep:  time series of discrete timesteps 
+%          edot: time series of energy generation 
+%          flux_end, the starting and ending points for each reaction flux.
+%          flux: time seris of the reaction fluxes in the network.
+%--------------------------------------------------------------------------
+
+  [~,aa,~,~,~,~] = read_nz_network(data_directory);
+  ny = length(aa);
+ 
+  sparse_filename = strcat(data_directory,'/sparse_ind');
+  file_id = fopen(sparse_filename,'rb');
+
+% Read Run Descriptions
+  record_length1=fread(file_id,1,'int32');
+  nnz    =fread(file_id,1,'int32');
+  record_length2=fread(file_id,1,'int32');
+  
+  record_length1=fread(file_id,1,'int32');
+  ridx   =fread(file_id,nnz,'int32');
+  cidx   =fread(file_id,nnz,'int32');
+  pb     =fread(file_id,ny+1,'int32');
+  record_length2=fread(file_id,1,'int32');
+  
+  record_length1=fread(file_id,1,'int32');
+  nrext  =fread(file_id,4,'int32');
+  record_length2=fread(file_id,1,'int32');
+  
+  ns1 = zeros(nrext(1),1);
+  ns2 = zeros(nrext(2),2);
+  ns3 = zeros(nrext(3),3);
+  ns4 = zeros(nrext(4),4);
+  
+  record_length1=fread(file_id,1,'int32');
+  ns1(:,1) =fread(file_id,nrext(1),'int32');
+  ns2(:,1) =fread(file_id,nrext(2),'int32');
+  ns2(:,2) =fread(file_id,nrext(2),'int32');
+  record_length2=fread(file_id,1,'int32');
+  
+  for i = 1:3
+      record_length1=fread(file_id,1,'int32');
+      ns3(:,i) =fread(file_id,nrext(3),'int32');
+      record_length2=fread(file_id,1,'int32');
+  end
+  
+  for i = 1:4
+      record_length1=fread(file_id,1,'int32');
+      ns4(:,i) =fread(file_id,nrext(4),'int32');
+      record_length2=fread(file_id,1,'int32');
+  end
+  
+  fclose(file_id);
+  
+end
+
diff --git a/tools/plot/matlab/view_flux_max.m b/tools/plot/matlab/view_flux_max.m
new file mode 100644
index 00000000..888d5a1a
--- /dev/null
+++ b/tools/plot/matlab/view_flux_max.m
@@ -0,0 +1,42 @@
+function [] = view_flux_max(filename,color_bins)
+%--------------------------------------------------------------------------
+%[] = view_flux_integrated(filename,color_bins)
+% Plots vectors for the time integrated reaction fluxes.  
+% Up to 6 sets of fluxes are plotted in different colors.
+% Inputs>  filename: file from which flux data is read.
+%          cycle_limit: The limits of integration, [1 0] inclues all.
+%          color_bins: array specifing the range of flux for each color. 
+% Outputs: None
+%--------------------------------------------------------------------------
+  font_size = 16;
+
+% Read TS file
+  [zz, aa, xmf, time, temperature, density, timestep, ~, flx_end, flx] = read_ts_file(filename);
+  nn=aa-zz;
+
+% Plot nuclear chart
+  [z_max,z_min,n_max,n_min,point_size] = draw_nz_background(zz,nn,font_size);
+
+% Identify starting and ending points for all fluxes
+  zt=zz(flx_end(:,1));  nt=nn(flx_end(:,1));
+  zp=zz(flx_end(:,2));  np=nn(flx_end(:,2));
+  
+% Find max fluxes
+  flx_write = zeros(size(flx,1),1);
+  [~,iwrite] = max(abs(flx),[],2);
+  for j = 1:size(flx,1)
+      flx_write(j) = flx(j,iwrite(j));
+  end
+  
+% Label Figure
+  cond_label={['Maximum flux ',num2str(1,'%5d'),' to ',num2str(size(time,2),'%5d')],...
+    ['Time =',num2str(1,'%8.3e'),' to ',num2str(size(time,2),'%8.3e'),' seconds']};
+   text(n_min+1, z_max-1.5,cond_label,'FontSize',14,'HorizontalAlignment','left');
+
+% Choose maximum flux
+  flx_max=max(abs(flx_write));
+
+% Draw flux vectors  
+  draw_nz_flux(flx_max,color_bins,flx_write,zt,nt,zp,np)
+  
+end
\ No newline at end of file