diff --git a/asgard_test.m b/asgard_test.m
index 16180898..a9ed5a90 100644
--- a/asgard_test.m
+++ b/asgard_test.m
@@ -265,3 +265,14 @@ function asgard_fokkerplanck1_5p1a_implicit_test(testCase)
 verifyLessThan(testCase,err,2.5e-2);
 end
 
+function asgard_BDF2_test(testCase)
+addpath(genpath(pwd));
+disp('Testing BDF2');
+[err,act_f,act_frs] = asgard(fokkerplanck1_4p1a,'lev',4,'deg',5,'timestep_method','BDF2', 'dt',0.01,'num_steps',1, 'quiet', true);
+[err1,act_f,act_frs] = asgard(fokkerplanck1_4p1a,'lev',4,'deg',5,'timestep_method','BDF2', 'dt',0.005,'num_steps',2,'quiet', true);
+[err2,act_f,act_frs] = asgard(fokkerplanck1_4p1a,'lev',4,'deg',5,'timestep_method','BDF2', 'dt',0.0025,'num_steps',4,'quiet', true);
+logslope = (log(err1) - log(err))/(log(1/0.1) - log(1/0.05));
+slopeerr = abs(logslope - 2.000);
+verifyLessThan(testCase,slopeerr,2.5e-2);
+end
+
diff --git a/runtime_defaults.m b/runtime_defaults.m
index 022a8ac5..b992809a 100644
--- a/runtime_defaults.m
+++ b/runtime_defaults.m
@@ -16,7 +16,7 @@
 default_use_oldhash = false;
 default_use_oldcoeffmat = false;
 default_timestep_method = 'RK3';
-valid_timestep_methods = {'BE','CN','ode15i','ode15s','ode45','RK3'};
+valid_timestep_methods = {'BE','CN','ode15i','ode15s','ode45','RK3','BDF2'};
 check_timestep_method = @(x) any(validatestring(x,valid_timestep_methods));
 default_time_independent_A = false;
 default_many_solution_capable = false;
@@ -125,7 +125,7 @@
 
 opts.timestep_method = input_parser.Results.timestep_method;
 opts.build_A = false;
-if sum(strcmp(opts.timestep_method,{'BE','CN'}))>0
+if sum(strcmp(opts.timestep_method,{'BE','CN','BDF2'}))>0
     opts.build_A = true;
 end
 
diff --git a/time_advance.m b/time_advance.m
index 07a9c15c..01ea9980 100644
--- a/time_advance.m
+++ b/time_advance.m
@@ -1,17 +1,20 @@
-function f = time_advance(pde,opts,A_data,f,t,dt,deg,hash_table,Vmax,Emax)
+function f1 = time_advance(pde,opts,A_data,f0,t,dt,deg,hash_table,Vmax,Emax)
 
 if strcmp(opts.timestep_method,'BE')
     % Backward Euler (BE) first order
-    f = backward_euler(pde,opts,A_data,f,t,dt,deg,hash_table,Vmax,Emax);
+    f1 = backward_euler(pde,opts,A_data,f0,t,dt,deg,hash_table,Vmax,Emax);
 elseif strcmp(opts.timestep_method,'CN')
     % Crank Nicolson (CN) second order
-    f = crank_nicolson(pde,opts,A_data,f,t,dt,deg,hash_table,Vmax,Emax);
+    f1 = crank_nicolson(pde,opts,A_data,f0,t,dt,deg,hash_table,Vmax,Emax);
+elseif strcmp(opts.timestep_method,'BDF2')
+    % Backward Difference Forumla (BDF) second order
+    f1 = bdf2(pde,opts,A_data,f0,t,dt,deg,hash_table,Vmax,Emax);
 elseif sum(strcmp(opts.timestep_method,{'ode15i','ode15s','ode45'}))>0
     % One of the atlab ODE integrators.
-    f = ODEm(pde,opts,A_data,f,t,dt,deg,hash_table,Vmax,Emax);  
+    f1 = ODEm(pde,opts,A_data,f0,t,dt,deg,hash_table,Vmax,Emax);  
 else
     % RK3 TVD
-    f = RungeKutta3(pde,opts,A_data,f,t,dt,deg,hash_table,Vmax,Emax);
+    f1 = RungeKutta3(pde,opts,A_data,f0,t,dt,deg,hash_table,Vmax,Emax);
 end
 
 end
@@ -49,7 +52,7 @@
 if strcmp(opts.timestep_method,'ode45')
     
     if(~opts.quiet);disp('Using ode45');end
-    options = odeset('RelTol',1e-3,'AbsTol',1e-6,'Stats','off');
+    options = odeset('RelTol',1e-3,'AbsTol',1e-6,'Stats','on');
     [tout,fout] = ode45(@explicit_ode,[t0 t0+dt],f0,options);
     
 elseif strcmp(opts.timestep_method,'ode15s')
@@ -72,7 +75,7 @@
     
     % call ode15s
     if(~opts.quiet);disp('Using ode15s');end
-    options = odeset('RelTol',1e-4,'AbsTol',1e-3,...
+    options = odeset('RelTol',1e-4,'AbsTol',1e-6,...
         'Stats','on','OutputFcn',@odetpbar,'Refine',20);%,'Jacobian', J2);%'JPattern',S);
     [tout,fout] = ode15s(@explicit_ode,[t0 t0+dt],f0,options);
     
@@ -237,3 +240,51 @@
 end
 
 end
+
+
+function f2 = bdf2(pde,opts,A_data,f0,t,dt,deg,hash_table,Vmax,Emax)
+
+s1 = source_vector(pde,opts,hash_table,t+dt/2);
+s2 = source_vector(pde,opts,hash_table,t+dt);
+
+bc1 = boundary_condition_vector(pde,opts,hash_table,t+dt/2);
+bc2 = boundary_condition_vector(pde,opts,hash_table,t+dt);
+
+% %%
+% Apply any non-identity LHS mass matrix coefficient
+
+applyLHS = ~isempty(pde.termsLHS);
+
+%Obtain f1 through same process as Backward Euler
+if t == 0
+%    opts.timestep_method = 'ode45';
+%    f1 = ODEm(pde,opts,A_data,f0,t,dt,deg,hash_table,Vmax,Emax);
+    fac = 200;
+    for n=1:fac
+        f0 = backward_euler(pde,opts,A_data,f0,t,dt/fac,deg,hash_table,Vmax,Emax);
+        %f0 = RungeKutta3(pde,opts,A_data,f0,t,dt/fac,deg,hash_table,Vmax,Emax);
+       t = t + dt/fac;
+    end
+    f1 = f0;
+    save('f1.mat', 'f1');
+else
+    load('f1.mat', 'f1'); %Need to update this at the end
+end
+
+if applyLHS %Now we only need BDF2
+    %Obtain updated A and ALHS matrices
+    [~,A2,ALHS2] = apply_A(pde,opts,A_data,f1,deg);
+    I = eye(numel(diag(A2)));
+    AA2 = I - 2/3*dt*(ALHS2 \ A2);
+    b2 = 4*f1/3 - f0/3 + 2/3*dt*(ALHS2 \ (bc2 + s2));
+else
+    [~,A2] = apply_A(pde,opts,A_data,f1,deg);
+    I = eye(numel(diag(A2)));
+    AA2 = I -2/3*dt*A2;
+    b2 = 4*f1/3 - f0/3 + 2/3*dt*(bc2 + s2);
+end
+f2 = AA2 \ b2;
+f1 = f2 ;
+save('f1.mat', 'f1');
+
+end
\ No newline at end of file
diff --git a/two_scale/two_scale_rel_2.mat b/two_scale/two_scale_rel_2.mat
index be266753..9422cf56 100644
Binary files a/two_scale/two_scale_rel_2.mat and b/two_scale/two_scale_rel_2.mat differ