diff --git a/analyze_matrix.m b/analyze_matrix.m
new file mode 100644
index 00000000..8053f235
--- /dev/null
+++ b/analyze_matrix.m
@@ -0,0 +1,49 @@
+function analyze_matrix(A)
+
+fprintf('\n\nbegin matrix analysis\n');
+
+% these tests disabled - they are slow, and none of the problem
+% matrices tested are symmetric or positive definite
+
+% 1: test for symmetry
+is_sym = isequal(A,A');
+
+if is_sym
+    fprintf("--symmetric--\n");
+else
+    fprintf("--nonsymmetric--\n");
+end
+
+
+[~,flag] = chol(A, 'lower');
+
+if flag == 0    
+    fprintf("--symmetric positive definite--\n");
+else
+    fprintf("--not symmetric positive definite--\n");
+end
+
+ 
+% 3: spectrum
+[eigenvectors, eigenvalue_mat] = eig(A);
+eigenvalues = diag(eigenvalue_mat);
+figure
+plot(eigenvalues,'o','MarkerSize',12);
+saveas(gcf,'spectrum.png');
+
+
+% 4: condition number of eigenvect mat
+fprintf('eigen condition number: %d\n', cond(eigenvectors));
+
+% 5: condition number of matrix
+fprintf('matrix condition number: %d\n', cond(A));
+
+% 6: sparsity %
+fprintf('matrix density: %f\n', nnz(A)/numel(A));
+
+% 7: sparsity pattern
+spy(A);
+
+fprintf('end matrix analysis\n\n');
+
+end
\ No newline at end of file
diff --git a/apply_A.m b/apply_A.m
index a3113277..79466637 100644
--- a/apply_A.m
+++ b/apply_A.m
@@ -34,6 +34,7 @@
 if opts.implicit
     A = zeros(totalDOF,totalDOF); % Only filled if implicit
 end
+ALHS = 0;
 if num_terms_LHS > 0
     ALHS = zeros(totalDOF,totalDOF); % Only filled if non-identity LHS mass matrix
 end
diff --git a/asgard.m b/asgard.m
index 69283e88..21c3a193 100644
--- a/asgard.m
+++ b/asgard.m
@@ -10,10 +10,8 @@
 runtime_defaults
 
 %% Reset any persistent variables
-if opts.time_independent_A
-    clear time_advance
-end
-
+clear time_advance
+figure(1000);
 %% Check PDE
 pde = check_pde(pde,opts);
 
diff --git a/iterative_solve.m b/iterative_solve.m
new file mode 100644
index 00000000..693e9254
--- /dev/null
+++ b/iterative_solve.m
@@ -0,0 +1,68 @@
+function [x, iter, relres, relres_hist] = iterative_solve(pde, opts, A, b, tol, max_iter, restart, graph_residual, precond)
+%iterative solve driver
+
+% set options
+if ~exist('restart','var')
+    restart = [];
+else
+   if strcmpi(opts.solve_choice, 'BICGSTAB')
+       disp('non-restarting method')
+   end
+end
+
+if ~exist('tol','var') || isempty(tol)
+    tol = 1e-6; % matlab default
+end
+
+if ~exist('max_iter', 'var') || isempty(max_iter)
+    max_iter = size(A, 1);
+end
+
+% set up preconditioner
+if strcmpi(opts.preconditioner, 'JACOBI') %point jacobi
+    M = diag(diag(A));
+elseif strcmpi(opts.preconditioner, 'BLK_JACOBI') % block jacobi
+   blk_n = pde.deg^size(pde.dimensions,1);
+   nblks = size(A,1) / blk_n; 
+   kron_map = kron(eye(nblks), ones(blk_n, blk_n));
+   blk_list = A(kron_map ~= 0);
+   blocks = reshape(blk_list,blk_n,blk_n,size(A,1)/blk_n);
+   M = num2cell(blocks, [1,2]);
+   M = reshape(M, 1, size(M, 3));
+   M = blkdiag(M{:});
+elseif strcmpi(opts.preconditioner, 'CUSTOM')
+    if ~exist('precond', 'var')
+        error('selected custom preconditioning, did not supply one');
+    end
+    M = precond; %custom preconditioner
+else
+    M = [];
+end
+
+% do solve
+tic;
+if strcmpi(opts.solve_choice, 'GMRES')
+    [x,flag, relres, iter, relres_hist] = gmres(A, b, restart, tol, max_iter, M);
+    if flag == 0; fprintf('GMRES converged w iters\n'); disp(iter);
+    else; fprintf('GMRES failed with flag %d\n', flag); end
+elseif strcmpi(opts.solve_choice, 'BICGSTAB')
+    [x,flag, relres, iter, relres_hist] = bicgstab(A, b, tol, max_iter, M);
+    if flag == 0; fprintf('BIGCGSTAB converged w iters\n'); disp(iter);
+    else; fprintf('BICGSTAB failed with flag %d\n', flag); end
+else
+    error('solver choice not recongized!');
+end
+disp(['solve duration: ', num2str(toc), 's']);
+
+
+% graph performance
+if graph_residual
+    figure(2000);
+    semilogy(1:size(relres_hist,1),relres_hist/norm(b),'-o');
+    xlabel('Iteration number');
+    ylabel('Relative residual');
+    
+end
+    
+end
+
diff --git a/runtime_defaults.m b/runtime_defaults.m
index 73ff890d..482e9192 100644
--- a/runtime_defaults.m
+++ b/runtime_defaults.m
@@ -25,6 +25,14 @@
 default_adapt_threshold = 1e-1;
 default_refinement_method = 1;
 default_adapt_initial_condition = false;
+default_solve_choice = 'DIRECT';
+valid_solve_choices = {'DIRECT', 'GMRES', 'BICGSTAB'};
+check_solve_choice = @(x) any(validatestring(x,valid_solve_choices));
+default_do_analysis = false;
+default_preconditioner_choice = 'JACOBI';
+valid_preconditioner_choices = {'JACOBI', 'BLK_JACOBI', 'CUSTOM', 'NONE'};
+check_preconditioner_choice = @(x) any(validatestring(x,valid_preconditioner_choices));
+
 
 addRequired(input_parser, 'pde', @isstruct);
 addParameter(input_parser,'lev',default_lev, @isnumeric);
@@ -46,6 +54,9 @@
 addOptional(input_parser,'adapt_threshold',default_adapt_threshold, @isnumeric);
 addOptional(input_parser,'refinement_method',default_refinement_method, @isnumeric);
 addOptional(input_parser,'adapt_initial_condition',default_adapt_initial_condition,@islogical);
+addOptional(input_parser,'analyze_matrix',default_do_analysis, @islogical);
+addOptional(input_parser,'solve_choice',default_solve_choice, check_solve_choice);
+addOptional(input_parser,'preconditioner',default_preconditioner_choice, check_preconditioner_choice);
 
 if numel(varargin) == 0 && ~exist('pde','var')
     
@@ -135,6 +146,13 @@
 opts.adapt_threshold = input_parser.Results.adapt_threshold;
 opts.refinement_method = input_parser.Results.refinement_method;
 opts.adapt_initial_condition = input_parser.Results.adapt_initial_condition;
+opts.analyze_matrix = input_parser.Results.analyze_matrix;
+opts.solve_choice = input_parser.Results.solve_choice;
+opts.preconditioner = input_parser.Results.preconditioner;
+
+if opts.analyze_matrix && ~opts.implicit
+    error('cannot analyze matrix if not build for implicit stepping');
+end
 
 if opts.adapt
     opts.use_oldhash = false;
diff --git a/time_advance.m b/time_advance.m
index c5864f0f..f2b3ab63 100644
--- a/time_advance.m
+++ b/time_advance.m
@@ -81,6 +81,10 @@
 %% Backward Euler (first order implicit time advance)
 
 function f1 = backward_euler(pde,opts,A_data,f0,t,dt,deg,hash_table,Vmax,Emax)
+persistent analysis_done;
+if isempty(analysis_done)
+    analysis_done = false;
+end
 
 s0 = source_vector(pde,opts,hash_table,t+dt);
 bc0 = boundary_condition_vector(pde,opts,hash_table,t+dt);
@@ -104,8 +108,20 @@
     b = f0 + dt*(s0 + bc0);
 end
 
+if opts.analyze_matrix && ~analysis_done
+    analyze_matrix(AA);
+    analysis_done = true;
+end
+
+
+if strcmpi(opts.solve_choice, 'DIRECT')
+    f1 = AA\b; % Solve at each timestep
+else
+    restart = 200;
+    graph_residual=true;
+    f1 = iterative_solve(pde, opts, AA, b, [], [], restart, graph_residual); 
+end
 
-f1 = AA\b; % Solve at each timestep
 
 end
 
@@ -113,6 +129,10 @@
 %% Crank Nicolson (second order implicit time advance)
 
 function f1 = crank_nicolson(pde,opts,A_data,f0,t,dt,deg,hash_table,Vmax,Emax)
+persistent analysis_done;
+if isempty(analysis_done)
+    analysis_done = false;
+end
 
 s0 = source_vector(pde,opts,hash_table,t);
 s1 = source_vector(pde,opts,hash_table,t+dt);
@@ -151,6 +171,11 @@
     
     if ~opts.quiet; disp(['    rcond(AA) : ', num2str(rcond(AA))]); end
     
+    if opts.analyze_matrix && ~analysis_done
+        analyze_matrix(AA);
+        analysis_done = true;
+    end
+    
     AA_inv = inv(AA);
     
 %     f1 = AA\b; % Solve at each timestep