diff --git a/asgard.m b/asgard.m
index 69283e88..bfed80b9 100644
--- a/asgard.m
+++ b/asgard.m
@@ -148,6 +148,13 @@
 fval_realspace = wavelet_to_realspace(pde,opts,Meval,fval,hash_table);
 fval_realspace_analytic = get_analytic_realspace_solution_D(pde,opts,coord,t);
 
+test_func = 1; %test function 
+
+%Taking the moment of test_func with respect to the numerical distribution
+%function
+
+test_moment = moment_integral(pde, fval_realspace, test_func); 
+
 err_wavelet = sqrt(mean((fval(:) - fval_analytic(:)).^2));
 err_realspace = sqrt(mean((fval_realspace(:) - fval_realspace_analytic(:)).^2));
 if ~opts.quiet
diff --git a/moment_integral.m b/moment_integral.m
new file mode 100644
index 00000000..89b8b9e8
--- /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); 
+
+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
+
+ww = ww.*(xmax - xmin);
+
+%[x, w] = lgwt(deg, 0, h);
+
+%points = [];
+%num = 2^Lev*deg;
+
+
+%for i = 0:2^Lev(i)-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/two_scale/two_scale_rel_2.mat b/two_scale/two_scale_rel_2.mat
deleted file mode 100644
index be266753..00000000
Binary files a/two_scale/two_scale_rel_2.mat and /dev/null differ