diff --git a/BothLoops.cpp b/BothLoops.cpp
new file mode 100644
index 0000000..5ff7100
--- /dev/null
+++ b/BothLoops.cpp
@@ -0,0 +1,93 @@
+#include <Rcpp.h>
+#include <math.h>       // power
+#include <cmath>
+
+// [[Rcpp::export]]
+Rcpp::NumericVector BothLoops(Rcpp::NumericVector& vProductivity, Rcpp::NumericVector& vGridCapital,
+                              Rcpp::NumericMatrix& mOutput, Rcpp::NumericMatrix& mTransition,
+                              Rcpp::NumericVector& bbeta_, Rcpp::IntegerVector& nGridCapital_, 
+                              Rcpp::IntegerVector& nGridProductivity_){
+  
+  const int nGridCapital = nGridCapital_[0];
+  const int nGridProductivity = nGridProductivity_[0];
+  const double bbeta = bbeta_[0];
+  
+  double valueProvisional, valueHighSoFar, consumption, capitalChoice;
+  
+  int nProductivity, nProductivityNextPeriod, nCapital, nCapitalNextPeriod, gridCapitalNextPeriod;
+  
+  double maxDifference = 10.0, diff, diffHighSoFar;
+  double tolerance = 0.0000001;
+  int iteration = 0;
+  
+  Rcpp::NumericMatrix mValueFunctionNew(nGridCapital,nGridProductivity);
+  Rcpp::NumericMatrix mValueFunction(nGridCapital,nGridProductivity);
+  Rcpp::NumericMatrix mPolicyFunction(nGridCapital,nGridProductivity);
+  Rcpp::NumericMatrix expectedValueFunction(nGridCapital,nGridProductivity);
+  
+  while(maxDifference > tolerance) {
+    
+    for (nProductivity = 0;nProductivity<nGridProductivity;++nProductivity){
+      for (nCapital = 0;nCapital<nGridCapital;++nCapital){
+        expectedValueFunction(nCapital,nProductivity) = 0.0;
+        for (nProductivityNextPeriod = 0;nProductivityNextPeriod<nGridProductivity;++nProductivityNextPeriod){
+          expectedValueFunction(nCapital,nProductivity) += mTransition(nProductivity,nProductivityNextPeriod)*mValueFunction(nCapital,nProductivityNextPeriod);
+        }
+      }
+    }
+    
+    for (nProductivity = 0;nProductivity<nGridProductivity;++nProductivity){
+      
+      // We start from previous choice (monotonicity of policy function)
+      gridCapitalNextPeriod = 0;
+      
+      for (nCapital = 0;nCapital<nGridCapital;++nCapital){
+        
+        valueHighSoFar = -100000.0;
+        capitalChoice  = vGridCapital[0];
+        
+        for (nCapitalNextPeriod = gridCapitalNextPeriod;nCapitalNextPeriod<nGridCapital;++nCapitalNextPeriod){
+          
+          consumption = mOutput(nCapital,nProductivity)-vGridCapital[nCapitalNextPeriod];
+          valueProvisional = (1-bbeta)*log(consumption)+bbeta*expectedValueFunction(nCapitalNextPeriod,nProductivity);
+          
+          if (valueProvisional>valueHighSoFar){
+            valueHighSoFar = valueProvisional;
+            capitalChoice = vGridCapital[nCapitalNextPeriod];
+            gridCapitalNextPeriod = nCapitalNextPeriod;
+          }
+          else{
+            break; // We break when we have achieved the max
+          }
+          mValueFunctionNew(nCapital,nProductivity) = valueHighSoFar;
+          mPolicyFunction(nCapital,nProductivity) = capitalChoice;
+        }
+      }
+    }
+    
+    diffHighSoFar = -100000.0;
+    for (nProductivity = 0;nProductivity<nGridProductivity;++nProductivity){
+      for (nCapital = 0;nCapital<nGridCapital;++nCapital){
+        diff = std::abs(mValueFunction(nCapital,nProductivity)-mValueFunctionNew(nCapital,nProductivity));
+        if (diff>diffHighSoFar){
+          diffHighSoFar = diff;
+        }
+        mValueFunction(nCapital,nProductivity) = mValueFunctionNew(nCapital,nProductivity);
+      }
+    }
+    
+    maxDifference = diffHighSoFar;
+    
+    iteration = iteration+1;
+    if (iteration % 10 ==0 || iteration ==1){
+      Rcpp::Rcout <<"Iteration = "<<iteration<<", Sup Diff = "<< maxDifference<<"\n";
+    }
+  }
+  
+  Rcpp::Rcout <<"Iteration = "<<iteration<<", Sup Diff = "<<maxDifference<<"\n";
+  Rcpp::Rcout <<" \n";
+  Rcpp::Rcout <<"My check = "<< mPolicyFunction(999,2)<<"\n";
+  Rcpp::Rcout <<" \n";
+  
+  return mPolicyFunction;
+}
diff --git a/RBC_Rcpp_2.R b/RBC_Rcpp_2.R
new file mode 100644
index 0000000..41c4ea2
--- /dev/null
+++ b/RBC_Rcpp_2.R
@@ -0,0 +1,57 @@
+## 0. Housekeeping
+
+rm(list=ls())
+
+Rcpp::sourceCpp('BothLoops.cpp')
+
+ptm <- proc.time()
+
+##  1. Calibration
+
+aalpha <- 1/3     # Elasticity of output w.r.t. capital
+bbeta  <- 0.95    # Discount factor
+
+# Productivity values
+vProductivity <- c(0.9792, 0.9896, 1.0000, 1.0106, 1.0212)
+
+# Transition matrix
+mTransition <- c(0.9727, 0.0273, 0.0000, 0.0000, 0.0000,
+                 0.0041, 0.9806, 0.0153, 0.0000, 0.0000,
+                 0.0000, 0.0082, 0.9837, 0.0082, 0.0000,
+                 0.0000, 0.0000, 0.0153, 0.9806, 0.0041,
+                 0.0000, 0.0000, 0.0000, 0.0273, 0.9727)
+
+mTransition <- matrix(mTransition,nrow=5,ncol=5)
+mTransition <- t(mTransition)
+
+## 2. Steady State
+
+capitalSteadyState <- (aalpha*bbeta)^(1/(1-aalpha))
+outputSteadyState <- capitalSteadyState^aalpha
+consumptionSteadyState <- outputSteadyState-capitalSteadyState
+
+cat(" Output = ", outputSteadyState,", Capital = ",capitalSteadyState, ", Consumption = ", consumptionSteadyState,"\n") 
+cat(" \n")
+
+# We generate the grid of capital
+vGridCapital <- seq(0.5*capitalSteadyState, 1.5*capitalSteadyState, by = 0.00001)
+
+nGridCapital <- length(vGridCapital)
+nGridProductivity <- length(vProductivity)
+
+# 3. Required matrices and vectors
+
+mOutput <- matrix(0,nGridCapital,nGridProductivity)
+
+## 4. We pre-build output for each point in the grid
+
+mOutput = vGridCapital^aalpha%*%t(vProductivity)
+
+## 5. Main iteration
+
+mPolicyFunction <- BothLoops(
+  vProductivity,vGridCapital, mOutput, mTransition, bbeta, nGridCapital, 
+  nGridProductivity
+)
+
+cat(" Time = ", proc.time() - ptm,"\n")