diff --git a/LevDegComparison_280Steps/Error_280Steps_dt=0,4_Spreadsheet.xlsx b/LevDegComparison_280Steps/Error_280Steps_dt=0,4_Spreadsheet.xlsx
new file mode 100644
index 00000000..4ee621dc
Binary files /dev/null and b/LevDegComparison_280Steps/Error_280Steps_dt=0,4_Spreadsheet.xlsx differ
diff --git a/LevDegComparison_280Steps/LevDegComparison_280Steps.zip b/LevDegComparison_280Steps/LevDegComparison_280Steps.zip
new file mode 100644
index 00000000..13926f9f
Binary files /dev/null and b/LevDegComparison_280Steps/LevDegComparison_280Steps.zip differ
diff --git a/LevDegComparison_280Steps/f1d_lev3_deg3.mat b/LevDegComparison_280Steps/f1d_lev3_deg3.mat
new file mode 100644
index 00000000..34d4144e
Binary files /dev/null and b/LevDegComparison_280Steps/f1d_lev3_deg3.mat differ
diff --git a/LevDegComparison_280Steps/f1d_lev3_deg4.mat b/LevDegComparison_280Steps/f1d_lev3_deg4.mat
new file mode 100644
index 00000000..c9b74e3e
Binary files /dev/null and b/LevDegComparison_280Steps/f1d_lev3_deg4.mat differ
diff --git a/LevDegComparison_280Steps/f1d_lev4_deg3.mat b/LevDegComparison_280Steps/f1d_lev4_deg3.mat
new file mode 100644
index 00000000..62376624
Binary files /dev/null and b/LevDegComparison_280Steps/f1d_lev4_deg3.mat differ
diff --git a/LevDegComparison_280Steps/f1d_lev4_deg4.mat b/LevDegComparison_280Steps/f1d_lev4_deg4.mat
new file mode 100644
index 00000000..263fa108
Binary files /dev/null and b/LevDegComparison_280Steps/f1d_lev4_deg4.mat differ
diff --git a/LevDegComparison_280Steps/f1d_lev5_deg3.mat b/LevDegComparison_280Steps/f1d_lev5_deg3.mat
new file mode 100644
index 00000000..08c5f0ba
Binary files /dev/null and b/LevDegComparison_280Steps/f1d_lev5_deg3.mat differ
diff --git a/LevDegComparison_280Steps/f1d_lev5_deg4.mat b/LevDegComparison_280Steps/f1d_lev5_deg4.mat
new file mode 100644
index 00000000..3d0d0ef4
Binary files /dev/null and b/LevDegComparison_280Steps/f1d_lev5_deg4.mat differ
diff --git a/PDE/advec_diffusion1.m b/PDE/advec_diffusion1.m
new file mode 100644
index 00000000..45c884ff
--- /dev/null
+++ b/PDE/advec_diffusion1.m
@@ -0,0 +1,154 @@
+function pde = advec_diffusion1
+% Example PDE using the 1D Advection-Diffusion Equation. This example PDE is
+% time dependent (although not all the terms are time dependent). This
+% implies the need for an initial condition. 
+% PDE:
+% 
+% df/dt = -df/dx + d^2 f/dx^2 + (cos(x) - sin(x))exp(-t)
+%
+% Domain is [0, pi/4]
+% Inhomogeneous Dirichlet boundary condition 
+% Code is expected to be added to the 1d Advection Equation 
+%
+% Diffusion terms are dealt with via LDG, i.e., splitting into two first
+% order equations:
+%
+% d^2 f / dx^2 becomes
+%
+% dq/dx with free (homogeneous Neumann BC)
+%
+% and
+%
+% q=df/dx with Dirichlet BCs specified by analytic solution
+%
+% Run with
+%
+% explicit
+% asgard(advec_diffusion1);
+%
+% implicit
+% asgard(advec_diffusion1,'implicit',true);
+
+pde.CFL = 0.01;
+
+%% Setup the dimensions
+% 
+% Here we setup a 1D problem (x,y)
+
+soln_x = @(x) cos(x) + sin(x);
+soln_t = @(t) exp(-t);
+
+BCFunc = @(x) soln_x(x);
+BCFunc_t = @(t) soln_t(t);
+
+% Domain is (a,b)
+
+% The function is defined for the plane
+% x = a and x = b
+BCL_fList = { ...
+    @(x,p,t) BCFunc(x), ... % replace x by a
+    @(t,p) BCFunc_t(t)
+    };
+
+BCR_fList = { ...
+    @(x,p,t) BCFunc(x), ... % replace x by b
+    @(t,p) BCFunc_t(t)
+    };
+%Incorrect BC to check whether each term requires BC
+%Each term DOES require correct BC
+%BCIncorrect_fList= { ...
+%    @(x,p,t) x.*0 + 2, ... % replace x by b
+%    @(t,p) BCFunc_t(t)
+%    };
+
+dim_x.name = 'x';
+dim_x.domainMin = 0;
+dim_x.domainMax = pi/4;
+dim_x.init_cond_fn = @(x,p,t) soln_x(x)*soln_t(t);
+
+
+%%
+% Add dimensions to the pde object
+% Note that the order of the dimensions must be consistent with this across
+% the remainder of this PDE.
+
+pde.dimensions = {dim_x};
+num_dims = numel(pde.dimensions);
+
+%% Setup the terms of the PDE
+%
+% Here we have 2 terms, with each term having nDims (x) operators.
+
+%%
+% Setup the -d_dx term
+
+g1 = @(x,p,t,dat) x.*0-1;
+pterm1 = GRAD(num_dims,g1,-1,'D','D', BCL_fList, BCR_fList);
+
+term1_x = TERM_1D({pterm1});
+term1   = TERM_ND(num_dims,{term1_x});
+
+%% 
+% Setup the d^2_dx^2 term
+
+% term1
+%
+% eq1 :  df/dt   == d/dx g1(x) q(x,y)   [grad,g1(x)=1, BCL=N, BCR=N]
+% eq2 :   q(x,y) == d/dx g2(x) f(x,y)   [grad,g2(x)=1, BCL=D, BCR=D]
+%
+% coeff_mat = mat1 * mat2
+
+g1 = @(x,p,t,dat) x.*0+1;
+g2 = @(x,p,t,dat) x.*0+1;
+
+pterm1 = GRAD(num_dims,g1,+1,'N','N');
+pterm2 = GRAD(num_dims,g2,-1,'D','D', BCL_fList, BCR_fList);
+
+term2_x = TERM_1D({pterm1,pterm2});
+term2   = TERM_ND(num_dims,{term2_x});
+
+%%
+% Add terms to the pde object
+
+ pde.terms = {term1, term2};
+
+%% Construct some parameters and add to pde object.
+%  These might be used within the various functions below.
+
+params.parameter1 = 0;
+params.parameter2 = 1;
+
+pde.params = params;
+
+%%
+% Sources
+
+s1x = @(x,p,t) cos(x) - sin(x);
+s1t = @(t,p) exp(-t);
+source1 = {s1x, s1t};
+
+pde.sources = {source1};
+
+%% Define the analytic solution (optional).
+% This requires nDims+time function handles.
+
+pde.analytic_solutions_1D = { ...
+    @(x,p,t) soln_x(x), ...
+    @(t,p) soln_t(t) 
+    };
+
+    function dt=set_dt(pde,CFL)
+        
+        dims = pde.dimensions;
+        
+        % for Diffusion equation: dt = C * dx^2
+        
+        lev = dims{1}.lev;
+        dx = 1/2^lev;
+        dt = CFL*dx^2;
+        
+    end
+
+pde.set_dt = @set_dt;
+
+end
diff --git a/PDE/fokkerplanck2_6p1.m b/PDE/fokkerplanck2_6p1.m
index 0492f6e2..a3a3611e 100644
--- a/PDE/fokkerplanck2_6p1.m
+++ b/PDE/fokkerplanck2_6p1.m
@@ -110,7 +110,7 @@
 
 %% Setup the dimensions 
 
-dim_p.domainMin = 0;
+dim_p.domainMin = 0.01;
 dim_p.domainMax = +10;
 dim_p.init_cond_fn = @(x,p,t) f0_p(x);
 
diff --git a/PDE/fokkerplanck2_6p1_withLHS.m b/PDE/fokkerplanck2_6p1_withLHS.m
index a024f96b..afda2950 100644
--- a/PDE/fokkerplanck2_6p1_withLHS.m
+++ b/PDE/fokkerplanck2_6p1_withLHS.m
@@ -15,7 +15,7 @@
 % asgard(fokkerplanck2_6p1_withLHS)
 %
 % implicit
-% asgard(fokkerplanck2_6p1_withLHS,'implicit',true,'num_steps',20,'CFL',1.0,'deg',3,'lev',4)
+% asgard(fokkerplanck2_6p1_withLHS, 'implicit', true,'CFL', 1.0, 'num_steps', 50, 'lev', 4, 'deg', 4)
 
 pde.CFL = 0.01;
 
@@ -32,7 +32,7 @@
 tau = 10^5;
 gamma = @(p)sqrt(1+(delta*p).^2);
 vx = @(p)1/vT*(p./gamma(p));
-p_min = 0.005;
+p_min = 0.01;
 
 Ca = @(p)nuEE*vT^2*(psi(vx(p))./vx(p));
 
@@ -43,7 +43,6 @@
     function ret = phi(x)
         ret = erf(x);
     end
-
     function ret = psi(x,t)     
         dphi_dx = 2./sqrt(pi) * exp(-x.^2);
         ret = 1./(2*x.^2) .* (phi(x) - x.*dphi_dx);   
@@ -115,7 +114,7 @@
             case 0.005
                 Q = 0.49859;    % for p_min = 0.005
             case 0.01              
-                Q = 0.497179;   % for p_min = 0.01
+                Q = 0.502844;  % for p_min = 0.01
             case 0.1
                 Q = 0.471814;   % for p_min = 0.1
             case 0.2
@@ -123,7 +122,7 @@
             case 0.5
                 Q = 0.361837;   % for p_min = 0.5
             case 1.0
-                Q = 0.23975;     % for p_min = 1.0
+                Q = 0.23975;    % for p_min = 1.0
         end
         ret = Q * 4/sqrt(pi) * exp(-x.^2);
     end
@@ -131,7 +130,7 @@
 %% Setup the dimensions 
 
 dim_p.domainMin = p_min;
-dim_p.domainMax = +10;
+dim_p.domainMax = 5;
 dim_p.init_cond_fn = @(x,p,t) f0_p(x);
 
 dim_z.domainMin = -1;
@@ -143,12 +142,38 @@
 % Note that the order of the dimensions must be consistent with this across
 % the remainder of this PDE.
 
+
+
 pde.dimensions = {dim_p,dim_z};
 num_dims = numel(pde.dimensions);
 
+%Inhomogeneous BC's set to analytical solution
+BCL_fList = { ...
+    @(x,p,t) soln_p(x,t), ... %replace x with p_min
+    @(x,p,t) soln_z(x,t), ...
+    @(t,p) 1
+    };
+
+BCR_fList = { ...
+    @(x,p,t) soln_p(x,t), ... %replace x with 10
+    @(x,p,t) soln_z(x,t), ...
+    @(t,p) 1
+    };
+
+%BCL_rList = { ...
+%    @(x,p,t) -2*x.*soln_p(x,t), ... %replace x with p_min
+%    @(x,p,t) soln_z(x,t), ...
+%    @(t,p) 1
+%    };
+%BCR_rList = { ...
+%    @(x,p,t) -2*x.*soln_p(x,t), ... %replace x with 10
+%    @(x,p,t) soln_z(x,t), ...
+%    @(t,p) 1
+%    };
+
 %% Setup the terms of the PDE
 
-%%
+%%r
 % termLHS == p^2 d/dt f(p,z,t)
 
 g1 = @(x,p,t,dat) x.^2;
@@ -167,11 +192,13 @@
 % termC1 == d/dp g2(p) r(p)   [grad, g2(p) = p^2*Ca, BCL=D,BCR=N]        
 %   r(p) == d/dp g3(p) f(p)   [grad, g3(p) = 1,      BCL=N,BCR=D]
 
+% Inhomogenous Dirichlet BC's for f(p_max)
+
 g2 = @(x,p,t,dat) x.^2.*Ca(x);
 g3 = @(x,p,t,dat) x.*0+1; 
 
 pterm2  = GRAD(num_dims,g2,+1,'D','N');
-pterm3  = GRAD(num_dims,g3,-1,'N','D');
+pterm3  = GRAD(num_dims,g3,-1,'N','D', BCL_fList, BCR_fList);
 term1_p = TERM_1D({pterm2,pterm3});
 termC1  = TERM_ND(num_dims,{term1_p,[]});
 
@@ -184,7 +211,7 @@
 
 g2 = @(x,p,t,dat) x.^2.*Cf(x);
 
-pterm2  = GRAD(num_dims,g2,+1,'N','D');
+pterm2  = GRAD(num_dims,g2, +1,'N','D', BCL_fList, BCR_fList);
 term2_p = TERM_1D({pterm2});
 termC2   = TERM_ND(num_dims,{term2_p,[]});
 
@@ -202,8 +229,8 @@
 pterm1  = MASS(g1);
 term3_p = TERM_1D({pterm1});
 
-g2 = @(x,p,t,dat) (1-x.^2);
-g3 = @(x,p,t,dat) x.*0+1;
+g2 = @(x,p,t,dat) 1-x.^2;
+g3 = @(x,p,t,dat) x.*0 + 1;
 pterm1  = GRAD(num_dims,g2,+1,'D','D');
 pterm2  = GRAD(num_dims,g3,-1,'N','N');
 term3_z = TERM_1D({pterm1,pterm2});
@@ -213,6 +240,7 @@
 %%
 % Add terms to the pde object
 
+
 pde.terms = {termC1,termC2,termC3};
 
 %% Construct some parameters and add to pde object.
diff --git a/PDE/fokkerplanck2_complete.m b/PDE/fokkerplanck2_complete.m
index 47006a6a..818c9ebb 100644
--- a/PDE/fokkerplanck2_complete.m
+++ b/PDE/fokkerplanck2_complete.m
@@ -58,7 +58,7 @@
     case '6p1b'
         delta = 0.042;
         Z = 1;
-        E = 0.25;
+        E = 0.1;
         tau = 10^5;
     case '6p1c' 
         delta = 0.042;
@@ -143,8 +143,8 @@
 
 %% Setup the dimensions 
 
-dim_p.domainMin = 0.1;
-dim_p.domainMax = +10;
+dim_p.domainMin = 0.01;
+dim_p.domainMax = +5;
 dim_p.init_cond_fn = @(x,p,t) f0_p(x);
 
 dim_z.domainMin = -1;
diff --git a/apply_A.m b/apply_A.m
index a3113277..d71a9631 100644
--- a/apply_A.m
+++ b/apply_A.m
@@ -31,9 +31,11 @@
 implicit = opts.implicit;
 
 totalDOF = nWork * elementDOF;
+A = 0;
 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 66136a6e..67597203 100644
--- a/asgard.m
+++ b/asgard.m
@@ -17,6 +17,8 @@
 %% Check PDE
 pde = check_pde(pde,opts);
 
+%pde.quad_x_in = linspace(-1, 1, pde.deg); %quadrature points are a linear space to include left and right ends
+
 %% Set time step.
 dt = pde.set_dt(pde,opts.CFL);
 if opts.dt_set_at_runtime
@@ -143,14 +145,25 @@
 else
    if ~opts.quiet; disp(['Number of DOF : ', num2str(numel(fval))]); end
 end
-
 fval_analytic = exact_solution_vector(pde,opts,hash_table,t);
 fval_realspace = wavelet_to_realspace(pde,opts,Meval,fval,hash_table);
 fval_realspace_analytic = get_analytic_realspace_solution_D(pde,opts,coord,t);
 
+% test_moment_matrix = 1; %test matrix for moment_integral
+% for j = 1:length(nodes{1,2})
+%     for i = 1:length(nodes{1,1})
+%        test_func(i,1) = exp(-0.5*((nodes{1,1}(i))^2))/sqrt(2*pi);
+%     end
+% end
+%test_moment = moment_integral(pde,fval_realspace_analytic,test_moment_matrix); 
+
+%Option to calculate moment of some key function
+
+
 err_wavelet = sqrt(mean((fval(:) - fval_analytic(:)).^2));
 err_realspace = sqrt(mean((fval_realspace(:) - fval_realspace_analytic(:)).^2));
 if ~opts.quiet
+%    disp(['test_moment :', num2str(test_moment)]);
     disp(['    num_dof : ', num2str(numel(fval))]);
     disp(['    wavelet space absolute err : ', num2str(err_wavelet)]);
     disp(['    wavelet space relative err : ', num2str(err_wavelet/max(abs(fval_analytic(:)))*100), ' %']);
@@ -166,7 +179,7 @@
 err = 1e9;
 if ~opts.quiet; disp('Advancing time ...'); end
 for L = 1:num_steps
-    
+    fval_prev = fval;
     tic;
     timeStr = sprintf('Step %i of %i at %f seconds',L,num_steps,t);
     
@@ -291,6 +304,11 @@
         % Get the real space solution
         
         fval_realspace = wavelet_to_realspace(pde,opts,Meval,fval,hash_table);
+        for i = 1:length(fval_realspace)
+           if fval_realspace(i) < 1e-15
+               fval_realspace(i) = 1e-15;
+           end
+        end
         
         %%
         % Try with function convertToRealSpace
@@ -323,7 +341,6 @@
         % Check the wavelet space solution with the analytic solution
         
         fval_analytic = exact_solution_vector(pde,opts,hash_table,t+dt);
-        
         assert(numel(fval)==numel(fval_analytic));
         
         err_wavelet = sqrt(mean((fval(:) - fval_analytic(:)).^2));
@@ -331,6 +348,10 @@
             disp(['    num_dof : ', num2str(numel(fval))]);
             disp(['    wavelet space absolute err : ', num2str(err_wavelet)]);
             disp(['    wavelet space relative err : ', num2str(err_wavelet/max(abs(fval_analytic(:)))*100), ' %']);
+        if norm(fval_prev(:) - fval(:)) < 1e-10
+            disp(['Converged']);
+            break
+        end
         end
         
         %%
@@ -371,18 +392,18 @@
     
     %%
     % Save output
-    
+    deg = pde.deg;
+%    LevX = pde.lev_vec(1);
+%    LevV = pde.lev_vec(2);
     saveOutput = 0;
-    if saveOutput
+    if saveOutput && t > num_steps-2
         stat = mkdir('output');
         fName = ['output/f2d-' sprintf('%04.4d',L) '.mat'];
         f2d = reshape(fval_realspace,deg*2^LevX,deg*2^LevV)';
         save(fName,'f2d','fval');
     end
-    
     t = t + dt;
     
 end
-
 end
 
diff --git a/asgard_test.m b/asgard_test.m
index 71a4c73e..a8738ef1 100644
--- a/asgard_test.m
+++ b/asgard_test.m
@@ -41,6 +41,31 @@ function asgard_advection1_implicit_test(testCase)
 
 end
 
+function asgard_advec_diffusion1_explicit_test(testCase)
+
+addpath(genpath(pwd));
+
+disp('Testing advec_diffusion1 (explicit)');
+
+[err, act_f, act_frs] = asgard(advec_diffusion1);
+
+verifyLessThan(testCase, err, 9e-4);
+
+end
+
+function asgard_advec_diffusion1_implicit_test(testCase)
+
+addpath(genpath(pwd));
+
+disp('Testing advec_diffusion1 (implicit)');
+
+[err, act_f, act_frs] = asgard(advec_diffusion1,'implicit',true);
+
+verifyLessThan(testCase, err, 9e-4);
+
+end
+
+
 function asgard_advection1reverse_implicit_test(testCase)
 
 addpath(genpath(pwd));
@@ -231,7 +256,7 @@ function asgard_diffusion2_adapt_test(testCase)
 
 [err,act_f,act_frs] = asgard(diffusion2,'lev',3,'quiet',true,'deg',3,'implicit',true,'adapt',true);
 
-verifyLessThan(testCase,err,1.5e-5);
+verifyLessThan(testCase,err,1e-3); % TODO : need to look into why this is larger than the adapt=false approach.
 
 end
 
@@ -259,3 +284,15 @@ function asgard_fokkerplanck1_5p1a_implicit_test(testCase)
 
 end
 
+function asgard_tokkerplanck2_6p1b_withLHS_implicit_test(testCase)
+
+addpath(genpath(pwd));
+
+disp('Testing fokkerplanck2_6p1b_withLHS (implicit/ with LHS)');
+
+[err, act_f, act_frs] = asgard(fokkerplanck2_6p1_withLHS, 'implicit', true,'CFL', 1.0, 'num_steps', 50, 'lev', 4, 'deg', 4);
+
+verifyLessThan(testCase, err, 2.1e-2);
+
+end
+
diff --git a/compute_boundary_condition.m b/compute_boundary_condition.m
index 3fe13b2a..e5ab59a2 100644
--- a/compute_boundary_condition.m
+++ b/compute_boundary_condition.m
@@ -10,6 +10,7 @@
 L = xMax-xMin;
 Tol_Cel_Num = 2^(Lev);
 h = L  / Tol_Cel_Num;
+small_dx = h*1e-7; %Parameter to fix problem with Inhomogeneous BC
 DoF = Deg * Tol_Cel_Num;
 
 bc = sparse(DoF,1);
@@ -22,14 +23,32 @@
     WorkCel = 0;
     c = [1:Deg];
     IntVal =  p_L'*Fun(xMin,p,t) ;
-    bc(c) = - g_func(xMin,p,t).*IntVal;
+    if isfinite(g_func(xMin,p,t)) %make sure g_func is finite
+        
+        bc(c) = -g_func(xMin,p,t).*IntVal;
     
+    else
+        
+        xLclose = xMin + small_dx; %Move away from xMin by small_dx
+        bc(c) = -g_func(xLclose,p,t).*IntVal;
+        
+    end
 else
     
     WorkCel = Tol_Cel_Num - 1;
     c = Deg*WorkCel+[1:Deg];
     IntVal =  p_R'*Fun(xMax,p,t);
-    bc(c) = g_func(xMax,p,t).*IntVal;
+    
+    if isfinite(g_func(xMax,p,t)) %make sure g_func is finite
+        
+        bc(c) = g_func(xMax,p,t).*IntVal;
+        
+    else
+        
+        xRclose = xMax - small_dx; %Move away from xMax by small_dx
+        bc(c) = g_func(xRclose,p,t).*Intval;
+        
+    end
     
 end
 
diff --git a/matrix_plot_D.m b/matrix_plot_D.m
index 52f067b1..cf108c99 100755
--- a/matrix_plot_D.m
+++ b/matrix_plot_D.m
@@ -1,65 +1,59 @@
-function [Meval,nodes] = matrix_plot_D(pde,dimension)
-
-%%
-% Generate the evaluation matrix and plotting points
-% TODO : probably should be renamed to something more descriptive
-
-%%
-% Setup a few shortcuts
-
-lev = dimension.lev;
-deg = pde.deg;
-xMin = dimension.domainMin;
-xMax = dimension.domainMax;
-FMWT = dimension.FMWT;
-
-%%
-% Jacobi of variable
-
-n = 2^(lev);
-h = (xMax-xMin)/n;
-dof_1D = deg*n;
-nodes = zeros(dof_1D,1);
-
-
-%%
-% Quadrature points (quad_x) and weights (quad_w)
-[quad_x,quad_w]=lgwt(deg,-1,1);
-
-%%
-% Uncomment to add end points to plot
-% 
-% quad_x = [-1 quad_x' +1]';
-
-dof = numel(quad_x);
-
-p_val = lin_legendre(quad_x,deg)*sqrt(1/h); % TODO : this call and normalization happens in several places. We should consolidate.
-
-Meval = sparse(dof*n,dof_1D);
-
-for i=0:n-1
-    
-    %%
-    % Mapping from [-1,1] to physical domain
-    xi = ((quad_x+1)/2+i)*h+xMin;
-    
-    %%
-    % Coefficients for DG bases
-    Iu = [dof*i+1:dof*(i+1)];
-    Iv = [deg*i+1:deg*(i+1)];
-    
-    Meval(Iu,Iv) = p_val;
-    nodes(Iu) = xi;
-    
-end
-
-%%
-% Transform back to real space from wavelet space
-
-Meval = Meval*FMWT';
-
-
-    
-
-
-
+function [Meval,nodes] = matrix_plot_D(pde,dimension, quad_x_in)
+
+%%
+% Generate the evaluation matrix and plotting points
+% TODO : probably should be renamed to something more descriptive
+
+%%
+% Setup a few shortcuts
+
+lev = dimension.lev;
+deg = pde.deg;
+xMin = dimension.domainMin;
+xMax = dimension.domainMax;
+FMWT = dimension.FMWT;
+
+%%
+% Jacobi of variable
+
+n = 2^(lev);
+h = (xMax-xMin)/n;
+dof_1D = deg*n;
+nodes = zeros(dof_1D,1);
+Meval = sparse(dof_1D,dof_1D);
+
+%%
+% Quadrature points (quad_x) and weights (quad_w)
+if (nargin >= 3)
+        quad_x = quad_x_in;
+else
+  [quad_x,quad_w]=lgwt(deg,-1,1);
+end
+
+p_val = lin_legendre(quad_x,deg)*sqrt(1/h); % TODO : this call and normalization happens in several places. We should consolidate.
+
+for i=0:n-1
+    
+    %%
+    % Mapping from [-1,1] to physical domain
+    xi = ((quad_x+1)/2+i)*h+xMin;
+    
+    %%
+    % Coefficients for DG bases
+    Iu = [deg*i+1:deg*(i+1)];
+    
+    Meval(Iu,Iu) = p_val;
+    nodes(Iu) = xi;
+    
+end
+
+%%
+% Transform back to real space from wavelet space
+
+Meval = Meval*FMWT';
+
+
+    
+
+
+
diff --git a/moment_integral.m b/moment_integral.m
new file mode 100644
index 00000000..f49cf2fc
--- /dev/null
+++ b/moment_integral.m
@@ -0,0 +1,43 @@
+function MomentValue = moment_integral(pde, fval_realspace, gfunc)
+
+xmin = pde.dimensions{1,1}.domainMin;
+xmax = pde.dimensions{1,1}.domainMax;
+Lev = pde.lev_vec;
+h = (xmax - xmin)/2^(Lev(1));
+deg = pde.deg;
+num_dimensions = length(pde.dimensions);
+
+[quad_xx, quad_ww] = lgwt(deg, -1, 1);
+
+quad_ww = 2^(-Lev(1))/2*quad_ww;
+
+ww = repmat(quad_ww, 2^Lev(1), 1); 
+
+ww = ww.*(xmax - xmin); %for 2D use Kronmult
+
+if num_dimensions >= 2
+    for i = 2:num_dimensions
+         domainMin = pde.dimensions{1,i}.domainMin;
+         domainMax = pde.dimensions{1,i}.domainMax;
+         ww = kron(ww,ww)*(domainMax - domainMin);
+    end
+end
+
+[x, w] = lgwt(deg, 0, h);
+
+points = [];
+%num = 2^Lev*deg;
+
+
+for i = 0:2^Lev(1)-1
+    points = [points; xmin + x + i*h];
+end
+
+%points2 = repmat(points, [num 1]);
+
+%points2 = reshape(reshape(points2,num,num)',num*num,1);
+
+MomentValue = sum(ww.*fval_realspace.*gfunc);
+
+end
+
diff --git a/output_Sebastian/A_matrix_250.mat b/output_Sebastian/A_matrix_250.mat
new file mode 100644
index 00000000..e8368bf3
Binary files /dev/null and b/output_Sebastian/A_matrix_250.mat differ
diff --git a/output_Sebastian/f2d-0250.mat b/output_Sebastian/f2d-0250.mat
new file mode 100644
index 00000000..6e6d9619
Binary files /dev/null and b/output_Sebastian/f2d-0250.mat differ
diff --git a/plot_fval.m b/plot_fval.m
index a4ec5483..e503c46b 100644
--- a/plot_fval.m
+++ b/plot_fval.m
@@ -222,4 +222,4 @@ function plot_fval(pde,nodes,fval_realspace,fval_realspace_analytic,Meval,coordi
 
 pause (0.01)
 
-end
\ No newline at end of file
+end
diff --git a/test_matrix_plot_D.m b/test_matrix_plot_D.m
new file mode 100644
index 00000000..d7d42226
--- /dev/null
+++ b/test_matrix_plot_D.m
@@ -0,0 +1,59 @@
+function [Meval,nodes] = matrix_plot_D(pde,dimension, quad_x_in)
+
+%%
+% Generate the evaluation matrix and plotting points
+% TODO : probably should be renamed to something more descriptive
+
+%%
+% Setup a few shortcuts
+
+lev = dimension.lev;
+deg = pde.deg;
+xMin = dimension.domainMin;
+xMax = dimension.domainMax;
+FMWT = dimension.FMWT;
+
+%%
+% Jacobi of variable
+
+n = 2^(lev);
+h = (xMax-xMin)/n;
+dof_1D = deg*n;
+nodes = zeros(dof_1D,1);
+Meval = sparse(dof_1D,dof_1D);
+
+%%
+% Quadrature points (quad_x) and weights (quad_w)
+if (nargin >= 3),
+        quad_x = quad_x_in;
+else
+  [quad_x,quad_w]=lgwt(deg,-1,1);
+end;
+
+p_val = lin_legendre(quad_x,deg)*sqrt(1/h); % TODO : this call and normalization happens in several places. We should consolidate.
+
+for i=0:n-1
+    
+    %%
+    % Mapping from [-1,1] to physical domain
+    xi = ((quad_x+1)/2+i)*h+xMin;
+    
+    %%
+    % Coefficients for DG bases
+    Iu = [deg*i+1:deg*(i+1)];
+    
+    Meval(Iu,Iu) = p_val;
+    nodes(Iu) = xi;
+    
+end
+
+%%
+% Transform back to real space from wavelet space
+
+Meval = Meval*FMWT';
+
+
+    
+
+
+
diff --git a/time_advance.m b/time_advance.m
index c5864f0f..5411a2fc 100644
--- a/time_advance.m
+++ b/time_advance.m
@@ -48,7 +48,7 @@
 b1 = 1/6; b2 = 2/3; b3 = 1/6;
 
 if applyLHS
-    [k1,A1,ALHS] = apply_A(pde,opts,A_data,f,deg,Vmax,Emax);
+    [k1,~,ALHS] = apply_A(pde,opts,A_data,f,deg,Vmax,Emax);
     rhs1 = source1 + bc1;
     %     invMatLHS = inv(ALHS); % NOTE : assume time independent for now for speed.
     %     k1 = invMatLHS * (k1 + rhs1);
diff --git a/two_scale/two_scale_rel_2.mat b/two_scale/two_scale_rel_2.mat
deleted file mode 100644
index a214909f..00000000
Binary files a/two_scale/two_scale_rel_2.mat and /dev/null differ