diff --git a/preprocessor/ExprNode.cc b/preprocessor/ExprNode.cc
index d9183d7119c77dd44e918ceadc36415b7528ad91..f8445b519e7f7cd12ee03f8702786f5cf22c4c0d 100644
--- a/preprocessor/ExprNode.cc
+++ b/preprocessor/ExprNode.cc
@@ -2916,8 +2916,7 @@ TrinaryOpNode::eval_opcode(double v1, TrinaryOpcode op_code, double v2, double v
   switch (op_code)
     {
     case oNormcdf:
-      cerr << "NORMCDF: eval not implemented" << endl;
-      exit(EXIT_FAILURE);
+      return (0.5*(1+erf((v1-v2)/v3/M_SQRT2)));
     }
   // Suppress GCC warning
   exit(EXIT_FAILURE);
@@ -2977,9 +2976,6 @@ void
 TrinaryOpNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
                            const temporary_terms_type &temporary_terms) const
 {
-  // TrinaryOpNode not implemented for C output
-  assert(!IS_C(output_type));
-
   // If current node is a temporary term
   temporary_terms_type::const_iterator it = temporary_terms.find(const_cast<TrinaryOpNode *>(this));
   if (it != temporary_terms.end())
@@ -2991,15 +2987,29 @@ TrinaryOpNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
   switch (op_code)
     {
     case oNormcdf:
-      output << "normcdf(";
+      if (IS_C(output_type))
+        {
+          // In C, there is no normcdf() primitive, so use erf()
+          output << "(0.5*(1+erf(((";
+          arg1->writeOutput(output, output_type, temporary_terms);
+          output << ")-(";
+          arg2->writeOutput(output, output_type, temporary_terms);
+          output << "))/(";
+          arg3->writeOutput(output, output_type, temporary_terms);
+          output << ")/M_SQRT2)))";
+        }
+      else
+        {
+          output << "normcdf(";
+          arg1->writeOutput(output, output_type, temporary_terms);
+          output << ",";
+          arg2->writeOutput(output, output_type, temporary_terms);
+          output << ",";
+          arg3->writeOutput(output, output_type, temporary_terms);
+          output << ")";
+        }
       break;
     }
-  arg1->writeOutput(output, output_type, temporary_terms);
-  output << ",";
-  arg2->writeOutput(output, output_type, temporary_terms);
-  output << ",";
-  arg3->writeOutput(output, output_type, temporary_terms);
-  output << ")";
 }
 
 void
diff --git a/tests/Makefile.am b/tests/Makefile.am
index d2d5b09d5fabf5a9b22001f5899a2fb850712902..c82f401a7bfbb9f2b1c3477c214a97e3882d430c 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -11,6 +11,7 @@ OCTAVE_MODS = \
 	osr_example.mod \
 	ramsey.mod \
 	ramst_initval_file.mod \
+	ramst_normcdf.mod \
 	example1_varexo_det.mod \
 	predetermined_variables.mod \
 	block_bytecode/fs2000_simk.mod \
diff --git a/tests/ramst_normcdf.mod b/tests/ramst_normcdf.mod
new file mode 100644
index 0000000000000000000000000000000000000000..89505144d4f270a1a1a1efb3f0a55946e0744281
--- /dev/null
+++ b/tests/ramst_normcdf.mod
@@ -0,0 +1,45 @@
+// Tests the normcdf() function, in the static M-file, and in a dynamic C file
+
+var c k t;
+varexo x;
+
+parameters alph gam delt bet aa;
+alph=0.5;
+gam=0.5;
+delt=0.02;
+bet=0.05;
+aa=0.5;
+
+
+model(use_dll);
+c + k - aa*x*k(-1)^alph - (1-delt)*k(-1);
+c^(-gam) - (1+bet)^(-1)*(aa*alph*x(+1)*k^(alph-1) + 1 - delt)*c(+1)^(-gam);
+t = normcdf(x, 2, 3);
+end;
+
+initval;
+x = 1;
+k = ((delt+bet)/(1.0*aa*alph))^(1/(alph-1));
+c = aa*k^alph-delt*k;
+t = 0;
+end;
+
+steady;
+
+check;
+
+shocks;
+var x;
+periods 1;
+values 1.2;
+end;
+
+simul(periods=20);
+
+if (abs(oo_.steady_state(3) - normcdf(1, 2, 3)) > 1e-10)
+   error('Test failed in static M-file')
+end
+
+if (abs(oo_.endo_simul(3, 2) - normcdf(1.2, 2, 3)) > 1e-10)
+   error('Test failed in dynamic C file')
+end