new file mode 100644
--- /dev/null
+++ b/main-Laplace.cpp
@@ -0,0 +1,113 @@
+//Solves 
+//
+//       nabla^2 u = f on  Omega
+//               u = g on dOmega
+//
+//Input data is read from interior.matrix and boundary.matrix
+//files. Each row in interior.matrix corresponds to an interior
+//point of Omega, and similarly for boundary.matrix. Although it's not
+//used, also need a normal vector at each point of the boundary as
+//given by normals.matrix
+
+#include <iostream>
+#include <fstream>
+#include <list>
+#include <cmath>
+
+#include <boost/shared_ptr.hpp>
+
+#include "include/linalg.hpp"
+#include "include/error.hpp"
+#include "include/rbf.hpp"
+#include "include/diff_op.hpp"
+#include "include/bvp.hpp"
+#include "include/interpolator.hpp"
+#include "include/utils.hpp"
+#include "include/func.hpp"
+
+#define RBF_TYPE inverse_multiquadric //Replace by template later.
+
+using namespace std;
+using namespace linalg;
+using namespace rbf;
+using namespace bvp;
+
+//Some arbitrary interior function for Poisson  equation (change to
+//"return 0*p(1);", for example, for Laplace equation).
+class intr_func : public realfunc{
+public:
+  double at(const point& p) const{return 0*(p(1)+p(2));};
+};
+
+//Boundary function for Poisson equation
+class bdry_func : public realfunc{
+public:
+  //sin(theta)^3
+  double at(const point& p) const{return 
+      std::pow(std::sin(std::atan2(p(1),p(2))),3);};
+};
+
+void output_result(const interpolator<RBF_TYPE>& u);
+
+int main(){
+  using linalg::vector;
+
+  try{
+    gsl_set_error_handler(&error_handling::errorHandler);
+    
+    //Define the domain from given data
+    shared_ptr<domain> Omega(new domain("data/interior.matrix",
+                                        "data/boundary.matrix",
+					"data/normals.matrix"));
+
+    //Define the interior and boundary operators
+    shared_ptr<Laplacian> L(new Laplacian);
+    shared_ptr<dirichlet_op> B(new dirichlet_op);
+
+    //rhs of BVP
+    intr_func f;
+    bdry_func g;
+    
+    //Define the Poisson boundary value problem.
+    shared_ptr<linear_BVP2> Poisson_BVP (new linear_BVP2(Omega, L, B, f, g));
+
+    //Define interpolator and solve the BVP.
+    interpolator<RBF_TYPE> u(Poisson_BVP);
+
+    //Write result to file
+    output_result(u);
+
+  }
+  catch(error& exc){
+    utils::show_exception(exc);
+  }
+
+  return 0;
+}
+
+void output_result(const interpolator<RBF_TYPE>& u)
+{
+  //Output result as a 41x41 meshgrid (to be plotted, say, with
+  //Octave; show_sol script included to work with this data).
+  matrix sol(41,41);    
+    
+  double r, th;
+  const double pi = atan(1)*4;
+  point p(2);
+  th = 0;
+  for(size_t i = 1; i <= 41; i++){
+    r = 0;
+    for(size_t j = 1; j <= 41; j++){
+      p(1) = r*cos(th); 
+      p(2) = r*sin(th);
+      sol(i,j) = u(p); //Evaluate interpolator at point p
+      r += 1.0/40;
+    }
+
+    th += 2.0*pi/40;
+  }
+
+    
+  ofstream ofs("data/sol.matrix");
+  ofs << sol;
+}