diff --git a/preprocessor b/preprocessor
index 78881c9244a087e23937acbfe0139014595d18aa..8bfcc59af7331621e534c8ff58eb5afeadb070fe 160000
--- a/preprocessor
+++ b/preprocessor
@@ -1 +1 @@
-Subproject commit 78881c9244a087e23937acbfe0139014595d18aa
+Subproject commit 8bfcc59af7331621e534c8ff58eb5afeadb070fe
diff --git a/tests/.gitignore b/tests/.gitignore
index d14f64f9274241b950d51a602027584183e61701..39a67982e2bb9f7066190782df394983aa9bc9d7 100644
--- a/tests/.gitignore
+++ b/tests/.gitignore
@@ -94,6 +94,9 @@ wsOct
 !/kronecker/test_kron.m
 !/load_octave_packages.m
 !/ls2003/data_ca1.m
+!/matrix_notation/extFunNoDerivsMatrix.m
+!/matrix_notation/extFunFirstDerivMatrix.m
+!/matrix_notation/extFunWithFirstAndSecondDerivsMatrix.m
 !/measurement_errors/data_ca1.m
 !/measurement_errors/fs2000_corr_me_ml_mcmc/fsdat_simul.m
 !/missing/simulate_data_with_missing_observations.m
diff --git a/tests/Makefile.am b/tests/Makefile.am
index b8333e1428c2c14aa1b27699713b7bee8b6e1a3f..474735e37923a3e4640d3af903af4b6c633b7bef 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -396,7 +396,18 @@ MODFILES = \
 	bgp/nk-1/nk.mod \
 	bgp/ramsey-1/ramsey.mod \
 	dynare-command-options/ramst.mod \
-	particle/local_state_space_iteration_k_test.mod
+	particle/local_state_space_iteration_k_test.mod \
+	matrix_notation/deterministic_matrix.mod \
+	matrix_notation/deterministic_scalar.mod \
+	matrix_notation/stochastic_matrix.mod \
+	matrix_notation/stochastic_scalar.mod \
+	matrix_notation/benchmark.mod \
+	matrix_notation/no_deriv_given_matrix_dll.mod \
+	matrix_notation/no_deriv_given_matrix.mod \
+	matrix_notation/first_deriv_given_matrix_dll.mod \
+	matrix_notation/first_deriv_given_matrix.mod \
+	matrix_notation/first_and_2nd_deriv_given_matrix_dll.mod \
+	matrix_notation/first_and_2nd_deriv_given_matrix.mod
 
 ECB_MODFILES = \
 	var-expectations/1/example.mod \
@@ -777,6 +788,23 @@ trend-component-and-var-models/tcm8.o.trs: trend-component-and-var-models/tcm6.o
 discretionary_policy/dennis_1_estim.m.trs: discretionary_policy/dennis_1.m.trs
 discretionary_policy/dennis_1_estim.o.trs: discretionary_policy/dennis_1.o.trs
 
+matrix_notation/deterministic_scalar.m.trs: matrix_notation/deterministic_matrix.m.trs
+matrix_notation/deterministic_scalar.o.trs: matrix_notation/deterministic_matrix.o.trs
+matrix_notation/stochastic_scalar.m.trs: matrix_notation/stochastic_matrix.m.trs
+matrix_notation/stochastic_scalar.o.trs: matrix_notation/stochastic_matrix.o.trs
+matrix_notation/first_deriv_given_matrix.m.trs: matrix_notation/benchmark.m.trs
+matrix_notation/first_deriv_given_matrix.o.trs: matrix_notation/benchmark.o.trs
+matrix_notation/first_deriv_given_matrix_dll.m.trs: matrix_notation/benchmark.m.trs
+matrix_notation/first_deriv_given_matrix_dll.o.trs: matrix_notation/benchmark.o.trs
+matrix_notation/first_and_2nd_deriv_given_matrix.m.trs: matrix_notation/benchmark.m.trs
+matrix_notation/first_and_2nd_deriv_given_matrix.o.trs: matrix_notation/benchmark.o.trs
+matrix_notation/first_and_2nd_deriv_given_matrix_dll.m.trs: matrix_notation/benchmark.m.trs
+matrix_notation/first_and_2nd_deriv_given_matrix_dll.o.trs: matrix_notation/benchmark.o.trs
+matrix_notation/no_deriv_given_matrix.m.trs: matrix_notation/benchmark.m.trs
+matrix_notation/no_deriv_given_matrix.o.trs: matrix_notation/benchmark.o.trs
+matrix_notation/no_deriv_given_matrix_dll.m.trs: matrix_notation/benchmark.m.trs
+matrix_notation/no_deriv_given_matrix_dll.o.trs: matrix_notation/benchmark.o.trs
+
 observation_trends_and_prefiltering/MCMC: m/observation_trends_and_prefiltering/MCMC o/observation_trends_and_prefiltering/MCMC
 m/observation_trends_and_prefiltering/MCMC: $(patsubst %.mod, %.m.trs, $(filter observation_trends_and_prefiltering/MCMC/%.mod, $(MODFILES)))
 o/observation_trends_and_prefiltering/MCMC: $(patsubst %.mod, %.o.trs, $(filter observation_trends_and_prefiltering/MCMC/%.mod, $(MODFILES)))
@@ -985,6 +1013,10 @@ estimation/univariate: m/estimation/univariate o/estimation/univariate
 m/estimation/univariate: $(patsubst %.mod, %.m.trs, $(filter estimation/univariate/%.mod, $(MODFILES)))
 o/estimation/univariate: $(patsubst %.mod, %.o.trs, $(filter estimation/univariate/%.mod, $(MODFILES)))
 
+matrix_notation: m/matrix_notation o/matrix_notation
+m/matrix_notation: $(patsubst %.mod, %.m.trs, $(filter matrix_notation/%.mod, $(MODFILES)))
+o/matrix_notation: $(patsubst %.mod, %.o.trs, $(filter matrix_notation/%.mod, $(MODFILES)))
+
 # ECB files
 M_ECB_TRS_FILES = $(patsubst %.mod, %.m.trs, $(ECB_MODFILES))
 
@@ -1168,7 +1200,10 @@ EXTRA_DIST = \
 	histval_initval_file/my_assert.m \
 	histval_initval_file/ramst_data.xls \
 	histval_initval_file/ramst_data.xlsx \
-	histval_initval_file/ramst_initval_file_data.m
+	histval_initval_file/ramst_initval_file_data.m \
+	matrix_notation/extFunFirstDerivMatrix.m \
+	matrix_notation/extFunNoDerivsMatrix.m \
+	matrix_notation/extFunWithFirstAndSecondDerivsMatrix.m
 
 if ENABLE_MATLAB
 check-local: check-matlab
diff --git a/tests/matrix_notation/benchmark.mod b/tests/matrix_notation/benchmark.mod
new file mode 100644
index 0000000000000000000000000000000000000000..c33ae1358ffa1d51f63fbbee12077a2fab7decc8
--- /dev/null
+++ b/tests/matrix_notation/benchmark.mod
@@ -0,0 +1,44 @@
+/* This file is used as a benchmark against with the other tests with external
+   functions are compared.
+   It is almost the same as example1.mod, except that the shocks process is
+   nonlinear (in order to have a non-zero Hessian of the
+   external function) */
+
+var y, c, k, h, a, b;
+varexo e, u;
+parameters beta, alpha, delta, theta, psi, rho, tau;
+
+alpha = 0.36;
+rho   = 0.95;
+tau   = 0.025;
+beta  = 0.99;
+delta = 0.025;
+psi   = 0;
+theta = 2.95;
+
+phi   = 0.1;
+
+model;
+c*theta*h^(1+psi)=(1-alpha)*y;
+k = beta*(((exp(b)*c)/(exp(b(+1))*c(+1)))
+          *(exp(b(+1))*alpha*y(+1)+(1-delta)*k));
+y = exp(a)*(k(-1)^alpha)*(h^(1-alpha));
+k = exp(b)*(y-c)+(1-delta)*k(-1) ;
+a = rho*a(-1)^2+5*tau*b(-1)^2 + e;
+b = 3*tau/2*a(-1)^2+7*rho/3*b(-1)^2 + u;
+end;
+
+initval;
+y = 1.08068253095672;
+c = 0.80359242014163;
+h = 0.29175631001732;
+k = 11.08360443260358;
+end;
+
+shocks;
+var e; stderr 0.009;
+var u; stderr 0.009;
+var e, u = phi*0.009*0.009;
+end;
+
+stoch_simul;
diff --git a/tests/matrix_notation/deterministic_matrix.mod b/tests/matrix_notation/deterministic_matrix.mod
new file mode 100644
index 0000000000000000000000000000000000000000..ff72d56c030f3ca8edff99e9bd1f88ac129d9658
--- /dev/null
+++ b/tests/matrix_notation/deterministic_matrix.mod
@@ -0,0 +1,38 @@
+// Do not forget to update deterministic_scalar.mod when this file is modified
+
+var(rows=3) X, Y, Z;
+var t;
+varexo(rows=3) E, U, V;
+parameters(rows=3, cols=3) P;
+
+P = [.1, .2, -.3; -.1, .5, .1; .2, -.3, .9];
+
+model;
+  X = P*X(-1)+E;
+  transpose(Y) = transpose(Y(-1))*transpose(P)+transpose(U);
+  Z = V;
+  t = Y[1](-1) + Y[2](+1);
+end;
+
+shocks;
+  var E;
+  periods 1 2;
+  values [0.5; 0.6; 0.7] [0.2; 0.3; 0.4];
+
+  var U[1];
+  periods 3;
+  values 0.8;
+
+  var U[2];
+  periods 3;
+  values 0.7;
+
+  var V[:];
+  periods 1;
+  values 0.1;
+end;
+
+perfect_foresight_setup(periods=50);
+perfect_foresight_solver;
+
+write_latex_dynamic_model;
diff --git a/tests/matrix_notation/deterministic_scalar.mod b/tests/matrix_notation/deterministic_scalar.mod
new file mode 100644
index 0000000000000000000000000000000000000000..3a893acc49c066e99fe029d63705d4481cf0b1f3
--- /dev/null
+++ b/tests/matrix_notation/deterministic_scalar.mod
@@ -0,0 +1,70 @@
+// Scalar translation of deterministic_matrix.mod
+// The declaration orders are the same as those of the expanded matrix version
+var t X_1 X_2 X_3 Y_1 Y_2 Y_3 Z_1 Z_2 Z_3;
+varexo E_1 E_2 E_3 U_1 U_2 U_3 V_1 V_2 V_3;
+parameters P_1_1 P_1_2 P_1_3 P_2_1 P_2_2 P_2_3 P_3_1 P_3_2 P_3_3;
+
+P_1_1 = .1;
+P_1_2 = .2;
+P_1_3 = -.3;
+P_2_1 = -.1;
+P_2_2 = .5;
+P_2_3 = .1;
+P_3_1 = .2;
+P_3_2 = -.3;
+P_3_3 = .9;
+
+model;
+  X_1 = P_1_1*X_1(-1) + P_1_2*X_2(-1) + P_1_3*X_3(-1) + E_1;
+  X_2 = P_2_1*X_1(-1) + P_2_2*X_2(-1) + P_2_3*X_3(-1) + E_2;
+  X_3 = P_3_1*X_1(-1) + P_3_2*X_2(-1) + P_3_3*X_3(-1) + E_3;
+  Y_1 = P_1_1*Y_1(-1) + P_1_2*Y_2(-1) + P_1_3*Y_3(-1) + U_1;
+  Y_2 = P_2_1*Y_1(-1) + P_2_2*Y_2(-1) + P_2_3*Y_3(-1) + U_2;
+  Y_3 = P_3_1*Y_1(-1) + P_3_2*Y_2(-1) + P_3_3*Y_3(-1) + U_3;
+  Z_1 = V_1;
+  Z_2 = V_2;
+  Z_3 = V_3;
+  t = Y_1(-1) + Y_2(+1);
+end;
+
+shocks;
+  var E_1;
+  periods 1 2;
+  values 0.5 0.2;
+
+  var E_2;
+  periods 1 2;
+  values 0.6 0.3;
+
+  var E_3;
+  periods 1 2;
+  values 0.7 0.4;
+
+  var U_1;
+  periods 3;
+  values 0.8;
+
+  var U_2;
+  periods 3;
+  values 0.7;
+
+  var V_1;
+  periods 1;
+  values 0.1;
+
+  var V_2;
+  periods 1;
+  values 0.1;
+
+  var V_3;
+  periods 1;
+  values 0.1;
+end;
+
+perfect_foresight_setup(periods=50);
+perfect_foresight_solver;
+
+L = load('deterministic_matrix_results.mat');
+if max(max(abs(L.oo_.endo_simul - oo_.endo_simul))) > 1e-12
+  error('Failure in matrix expansion')
+end
diff --git a/tests/matrix_notation/extFunFirstDerivMatrix.m b/tests/matrix_notation/extFunFirstDerivMatrix.m
new file mode 100644
index 0000000000000000000000000000000000000000..68c27d3a883cc179615b364d8bcee17405ba7b04
--- /dev/null
+++ b/tests/matrix_notation/extFunFirstDerivMatrix.m
@@ -0,0 +1,4 @@
+function dy=extFunFirstDerivMatrix(a,b,c)
+dy = [ c(1)^2 0        5*c(2)^2 0        2*a(1)*c(1) 10*b(1)*c(2);
+       0      3*c(1)^2 0        7*c(2)^2 6*a(2)*c(1) 14*b(2)*c(2) ];
+end
diff --git a/tests/matrix_notation/extFunNoDerivsMatrix.m b/tests/matrix_notation/extFunNoDerivsMatrix.m
new file mode 100644
index 0000000000000000000000000000000000000000..e7541cee70badbb3424555eee5c1a45608049e1b
--- /dev/null
+++ b/tests/matrix_notation/extFunNoDerivsMatrix.m
@@ -0,0 +1,3 @@
+function y=extFunNoDerivsMatrix(a,b,c)
+y=[a(1), 5*b(1); 3*a(2) 7*b(2)]*c.^2;
+end
diff --git a/tests/matrix_notation/extFunWithFirstAndSecondDerivsMatrix.m b/tests/matrix_notation/extFunWithFirstAndSecondDerivsMatrix.m
new file mode 100644
index 0000000000000000000000000000000000000000..e7b227909c2ef3dfa7461db60fff6e0e3295c6f8
--- /dev/null
+++ b/tests/matrix_notation/extFunWithFirstAndSecondDerivsMatrix.m
@@ -0,0 +1,20 @@
+function [y dy d2y]=extFunWithFirstAndSecondDerivsMatrix(a,b,c)
+y=[a(1), 5*b(1); 3*a(2), 7*b(2)]*c.^2;
+
+dy = [ c(1)^2 0        5*c(2)^2 0        2*a(1)*c(1) 10*b(1)*c(2);
+       0      3*c(1)^2 0        7*c(2)^2 6*a(2)*c(1) 14*b(2)*c(2) ];
+
+d2y=zeros(2,6,6);
+d2y(1,1,5) = 2*c(1);
+d2y(1,3,6) = 10*c(2);
+d2y(1,5,1) = 2*c(1);
+d2y(1,5,5) = 2*a(1);
+d2y(1,6,3) = 10*c(2);
+d2y(1,6,6) = 10*b(1);
+d2y(2,2,5) = 6*c(1);
+d2y(2,4,6) = 14*c(2);
+d2y(2,5,2) = 6*c(1);
+d2y(2,5,5) = 6*a(2);
+d2y(2,6,4) = 14*c(2);
+d2y(2,6,6) = 14*b(2);
+end
diff --git a/tests/matrix_notation/first_and_2nd_deriv_given_matrix.mod b/tests/matrix_notation/first_and_2nd_deriv_given_matrix.mod
new file mode 100644
index 0000000000000000000000000000000000000000..120ef1043a2c0ec03f482f72ab022816acb092bb
--- /dev/null
+++ b/tests/matrix_notation/first_and_2nd_deriv_given_matrix.mod
@@ -0,0 +1,71 @@
+var y, c, k, h;
+var(rows = 2) A; // A == [a; b]
+varexo(rows = 2) E; // E == [e; u]
+parameters beta, alpha, delta, theta, psi;
+parameters(rows=2,cols=2) P;
+
+alpha = 0.36;
+rho   = 0.95;
+tau   = 0.025;
+beta  = 0.99;
+delta = 0.025;
+psi   = 0;
+theta = 2.95;
+
+phi   = 0.1;
+
+P = [rho, tau; tau/2, rho/3];
+
+external_function(nargs=3, name=extFunWithFirstAndSecondDerivsMatrix, in_arg_dim=[2,2,2],
+                  out_arg_dim=2, first_deriv_provided, second_deriv_provided);
+
+model;
+[endogenous='c',name='law of motion of capital']
+c*theta*h^(1+psi)=(1-alpha)*y;
+k = beta*(((exp(A[2])*c)/(exp(A[2](+1))*c(+1)))
+          *(exp(A[2](+1))*alpha*y(+1)+(1-delta)*k));
+y = exp(A[1])*(k(-1)^alpha)*(h^(1-alpha));
+k = exp(A[2])*(y-c)+(1-delta)*k(-1) ;
+A = extFunWithFirstAndSecondDerivsMatrix(P[:,1], transpose(P[:,2]), A(-1)) + E;
+end;
+
+initval;
+y = 1.08068253095672;
+c = 0.80359242014163;
+h = 0.29175631001732;
+k = 11.08360443260358;
+
+A[1] = 0;
+A[2] = 0;
+E[:] = 0;
+end;
+
+shocks;
+var E[:]; stderr 0.009;
+var E[1], E[2] = phi*0.009*0.009;
+end;
+
+stoch_simul;
+
+L = load('benchmark_results.mat');
+if max(max(abs(L.oo_.dr.ghu - oo_.dr.ghu))) > 1e-12
+  error('Failure in matrix notation with external function')
+end
+if max(max(abs(L.oo_.dr.ghx - oo_.dr.ghx))) > 1e-12
+  error('Failure in matrix notation with external function')
+end
+if max(max(abs(L.oo_.dr.ghxu - oo_.dr.ghxu))) > 1e-12
+  error('Failure in matrix notation with external function')
+end
+if max(max(abs(L.oo_.dr.ghxx - oo_.dr.ghxx))) > 1e-12
+  error('Failure in matrix notation with external function')
+end
+if max(max(abs(L.oo_.dr.ghuu - oo_.dr.ghuu))) > 1e-12
+  error('Failure in matrix notation with external function')
+end
+if max(max(abs(L.oo_.dr.ghs2 - oo_.dr.ghs2))) > 1e-12
+  error('Failure in matrix notation with external function')
+end
+if max(max(abs(L.oo_.var - oo_.var))) > 1e-12
+  error('Failure in matrix notation with external function')
+end
diff --git a/tests/matrix_notation/first_and_2nd_deriv_given_matrix_dll.mod b/tests/matrix_notation/first_and_2nd_deriv_given_matrix_dll.mod
new file mode 100644
index 0000000000000000000000000000000000000000..570cf021dedddde43c01f11f0d8e5e3761a79f98
--- /dev/null
+++ b/tests/matrix_notation/first_and_2nd_deriv_given_matrix_dll.mod
@@ -0,0 +1,71 @@
+var y, c, k, h;
+var(rows = 2) A; // A == [a; b]
+varexo(rows = 2) E; // E == [e; u]
+parameters beta, alpha, delta, theta, psi;
+parameters(rows=2,cols=2) P;
+
+alpha = 0.36;
+rho   = 0.95;
+tau   = 0.025;
+beta  = 0.99;
+delta = 0.025;
+psi   = 0;
+theta = 2.95;
+
+phi   = 0.1;
+
+P = [rho, tau; tau/2, rho/3];
+
+external_function(nargs=3, name=extFunWithFirstAndSecondDerivsMatrix, in_arg_dim=[2,2,2],
+                  out_arg_dim=2, first_deriv_provided, second_deriv_provided);
+
+model(use_dll);
+[endogenous='c',name='law of motion of capital']
+c*theta*h^(1+psi)=(1-alpha)*y;
+k = beta*(((exp(A[2])*c)/(exp(A[2](+1))*c(+1)))
+          *(exp(A[2](+1))*alpha*y(+1)+(1-delta)*k));
+y = exp(A[1])*(k(-1)^alpha)*(h^(1-alpha));
+k = exp(A[2])*(y-c)+(1-delta)*k(-1) ;
+A = extFunWithFirstAndSecondDerivsMatrix(P[:,1], transpose(P[:,2]), A(-1)) + E;
+end;
+
+initval;
+y = 1.08068253095672;
+c = 0.80359242014163;
+h = 0.29175631001732;
+k = 11.08360443260358;
+
+A[1] = 0;
+A[2] = 0;
+E[:] = 0;
+end;
+
+shocks;
+var E[:]; stderr 0.009;
+var E[1], E[2] = phi*0.009*0.009;
+end;
+
+stoch_simul;
+
+L = load('benchmark_results.mat');
+if max(max(abs(L.oo_.dr.ghu - oo_.dr.ghu))) > 1e-12
+  error('Failure in matrix notation with external function')
+end
+if max(max(abs(L.oo_.dr.ghx - oo_.dr.ghx))) > 1e-12
+  error('Failure in matrix notation with external function')
+end
+if max(max(abs(L.oo_.dr.ghxu - oo_.dr.ghxu))) > 1e-12
+  error('Failure in matrix notation with external function')
+end
+if max(max(abs(L.oo_.dr.ghxx - oo_.dr.ghxx))) > 1e-12
+  error('Failure in matrix notation with external function')
+end
+if max(max(abs(L.oo_.dr.ghuu - oo_.dr.ghuu))) > 1e-12
+  error('Failure in matrix notation with external function')
+end
+if max(max(abs(L.oo_.dr.ghs2 - oo_.dr.ghs2))) > 1e-12
+  error('Failure in matrix notation with external function')
+end
+if max(max(abs(L.oo_.var - oo_.var))) > 1e-12
+  error('Failure in matrix notation with external function')
+end
diff --git a/tests/matrix_notation/first_deriv_given_matrix.mod b/tests/matrix_notation/first_deriv_given_matrix.mod
new file mode 100644
index 0000000000000000000000000000000000000000..431df3d64395ef2c97e1a4a6209a8c401eb7c136
--- /dev/null
+++ b/tests/matrix_notation/first_deriv_given_matrix.mod
@@ -0,0 +1,72 @@
+// --+ options: json=compute +--
+
+var y, c, k, h;
+var(rows = 2) A; // A == [a; b]
+varexo(rows = 2) E; // E == [e; u]
+parameters beta, alpha, delta, theta, psi;
+parameters(rows=2,cols=2) P;
+
+alpha = 0.36;
+rho   = 0.95;
+tau   = 0.025;
+beta  = 0.99;
+delta = 0.025;
+psi   = 0;
+theta = 2.95;
+
+phi   = 0.1;
+
+P = [rho, tau; tau/2, rho/3];
+
+external_function(nargs=3, name=extFunNoDerivsMatrix, in_arg_dim=[2,2,2],
+                  out_arg_dim=2, first_deriv_provided=extFunFirstDerivMatrix);
+
+model;
+c*theta*h^(1+psi)=(1-alpha)*y;
+k = beta*(((exp(A[2])*c)/(exp(A[2](+1))*c(+1)))
+          *(exp(A[2](+1))*alpha*y(+1)+(1-delta)*k));
+y = exp(A[1])*(k(-1)^alpha)*(h^(1-alpha));
+k = exp(A[2])*(y-c)+(1-delta)*k(-1) ;
+A = extFunNoDerivsMatrix(P[:,1], transpose(P[:,2]), A(-1)) + E;
+end;
+
+initval;
+y = 1.08068253095672;
+c = 0.80359242014163;
+h = 0.29175631001732;
+k = 11.08360443260358;
+
+A[1] = 0;
+A[2] = 0;
+E[:] = 0;
+end;
+
+shocks;
+var E[:]; stderr 0.009;
+var E[1], E[2] = phi*0.009*0.009;
+end;
+
+stoch_simul;
+
+L = load('benchmark_results.mat');
+if max(max(abs(L.oo_.dr.ghu - oo_.dr.ghu))) > 1e-12
+  error('Failure in matrix notation with external function')
+end
+if max(max(abs(L.oo_.dr.ghx - oo_.dr.ghx))) > 1e-12
+  error('Failure in matrix notation with external function')
+end
+if max(max(abs(L.oo_.dr.ghxu - oo_.dr.ghxu))) > 1e-12
+  error('Failure in matrix notation with external function')
+end
+if max(max(abs(L.oo_.dr.ghxx - oo_.dr.ghxx))) > 1e-12
+  error('Failure in matrix notation with external function')
+end
+if max(max(abs(L.oo_.dr.ghuu - oo_.dr.ghuu))) > 1e-12
+  error('Failure in matrix notation with external function')
+end
+if max(max(abs(L.oo_.dr.ghs2 - oo_.dr.ghs2))) > 1e-12
+  error('Failure in matrix notation with external function')
+end
+if max(max(abs(L.oo_.var - oo_.var))) > 1e-12
+  error('Failure in matrix notation with external function')
+end
diff --git a/tests/matrix_notation/first_deriv_given_matrix_dll.mod b/tests/matrix_notation/first_deriv_given_matrix_dll.mod
new file mode 100644
index 0000000000000000000000000000000000000000..070195c210c4d1c4db8a3b3a6de874f4189050e9
--- /dev/null
+++ b/tests/matrix_notation/first_deriv_given_matrix_dll.mod
@@ -0,0 +1,73 @@
+// --+ options: json=compute +--
+
+var y, c, k, h;
+var(rows = 2) A; // A == [a; b]
+varexo(rows = 2) E; // E == [e; u]
+parameters beta, alpha, delta, theta, psi;
+parameters(rows=2,cols=2) P;
+
+alpha = 0.36;
+rho   = 0.95;
+tau   = 0.025;
+beta  = 0.99;
+delta = 0.025;
+psi   = 0;
+theta = 2.95;
+
+phi   = 0.1;
+
+P = [rho, tau; tau/2, rho/3];
+
+external_function(nargs=3, name=extFunNoDerivsMatrix, in_arg_dim=[2,2,2],
+                  out_arg_dim=2, first_deriv_provided=extFunFirstDerivMatrix);
+
+model(use_dll);
+c*theta*h^(1+psi)=(1-alpha)*y;
+k = beta*(((exp(A[2])*c)/(exp(A[2](+1))*c(+1)))
+          *(exp(A[2](+1))*alpha*y(+1)+(1-delta)*k));
+y = exp(A[1])*(k(-1)^alpha)*(h^(1-alpha));
+k = exp(A[2])*(y-c)+(1-delta)*k(-1) ;
+A = extFunNoDerivsMatrix(P[:,1], transpose(P[:,2]), A(-1)) + E;
+end;
+
+initval;
+y = 1.08068253095672;
+c = 0.80359242014163;
+h = 0.29175631001732;
+k = 11.08360443260358;
+
+A[1] = 0;
+A[2] = 0;
+E[:] = 0;
+
+end;
+
+shocks;
+var E[:]; stderr 0.009;
+var E[1], E[2] = phi*0.009*0.009;
+end;
+
+stoch_simul;
+
+L = load('benchmark_results.mat');
+if max(max(abs(L.oo_.dr.ghu - oo_.dr.ghu))) > 1e-12
+  error('Failure in matrix notation with external function')
+end
+if max(max(abs(L.oo_.dr.ghx - oo_.dr.ghx))) > 1e-12
+  error('Failure in matrix notation with external function')
+end
+if max(max(abs(L.oo_.dr.ghxu - oo_.dr.ghxu))) > 1e-12
+  error('Failure in matrix notation with external function')
+end
+if max(max(abs(L.oo_.dr.ghxx - oo_.dr.ghxx))) > 1e-12
+  error('Failure in matrix notation with external function')
+end
+if max(max(abs(L.oo_.dr.ghuu - oo_.dr.ghuu))) > 1e-12
+  error('Failure in matrix notation with external function')
+end
+if max(max(abs(L.oo_.dr.ghs2 - oo_.dr.ghs2))) > 1e-12
+  error('Failure in matrix notation with external function')
+end
+if max(max(abs(L.oo_.var - oo_.var))) > 1e-12
+  error('Failure in matrix notation with external function')
+end
diff --git a/tests/matrix_notation/no_deriv_given_matrix.mod b/tests/matrix_notation/no_deriv_given_matrix.mod
new file mode 100644
index 0000000000000000000000000000000000000000..08c22198ffa5c313da9c1f24db60884ad4ee2ed0
--- /dev/null
+++ b/tests/matrix_notation/no_deriv_given_matrix.mod
@@ -0,0 +1,71 @@
+// --+ options: json=compute +--
+
+var y, c, k, h;
+var(rows = 2) A; // A == [a; b]
+varexo(rows = 2) E; // E == [e; u]
+parameters beta, alpha, delta, theta, psi;
+parameters(rows=2,cols=2) P;
+
+alpha = 0.36;
+rho   = 0.95;
+tau   = 0.025;
+beta  = 0.99;
+delta = 0.025;
+psi   = 0;
+theta = 2.95;
+
+phi   = 0.1;
+
+P = [rho, tau; tau/2, rho/3];
+
+external_function(nargs=3, name=extFunNoDerivsMatrix, in_arg_dim=[2,2,2], out_arg_dim=2);
+
+model;
+c*theta*h^(1+psi)=(1-alpha)*y;
+k = beta*(((exp(A[2])*c)/(exp(A[2](+1))*c(+1)))
+          *(exp(A[2](+1))*alpha*y(+1)+(1-delta)*k));
+y = exp(A[1])*(k(-1)^alpha)*(h^(1-alpha));
+k = exp(A[2])*(y-c)+(1-delta)*k(-1) ;
+A = extFunNoDerivsMatrix(P[:,1], transpose(P[:,2]), A(-1)) + E;
+end;
+
+initval;
+y = 1.08068253095672;
+c = 0.80359242014163;
+h = 0.29175631001732;
+k = 11.08360443260358;
+
+A[1] = 0;
+A[2] = 0;
+E[:] = 0;
+end;
+
+shocks;
+var E[:]; stderr 0.009;
+var E[1], E[2] = phi*0.009*0.009;
+end;
+
+stoch_simul;
+
+L = load('benchmark_results.mat');
+if max(max(abs(L.oo_.dr.ghu - oo_.dr.ghu))) > 1e-12
+  error('Failure in matrix notation with external function')
+end
+if max(max(abs(L.oo_.dr.ghx - oo_.dr.ghx))) > 1e-12
+  error('Failure in matrix notation with external function')
+end
+if max(max(abs(L.oo_.dr.ghxu - oo_.dr.ghxu))) > 1e-12
+  error('Failure in matrix notation with external function')
+end
+if max(max(abs(L.oo_.dr.ghxx - oo_.dr.ghxx))) > 1e-12
+  error('Failure in matrix notation with external function')
+end
+if max(max(abs(L.oo_.dr.ghuu - oo_.dr.ghuu))) > 1e-12
+  error('Failure in matrix notation with external function')
+end
+if max(max(abs(L.oo_.dr.ghs2 - oo_.dr.ghs2))) > 1e-12
+  error('Failure in matrix notation with external function')
+end
+if max(max(abs(L.oo_.var - oo_.var))) > 1e-12
+  error('Failure in matrix notation with external function')
+end
diff --git a/tests/matrix_notation/no_deriv_given_matrix_dll.mod b/tests/matrix_notation/no_deriv_given_matrix_dll.mod
new file mode 100644
index 0000000000000000000000000000000000000000..ccbedff1c778bf66748e7e9c5ec1074ea190fab3
--- /dev/null
+++ b/tests/matrix_notation/no_deriv_given_matrix_dll.mod
@@ -0,0 +1,71 @@
+// --+ options: json=compute +--
+
+var y, c, k, h;
+var(rows = 2) A; // A == [a; b]
+varexo(rows = 2) E; // E == [e; u]
+parameters beta, alpha, delta, theta, psi;
+parameters(rows=2,cols=2) P;
+
+alpha = 0.36;
+rho   = 0.95;
+tau   = 0.025;
+beta  = 0.99;
+delta = 0.025;
+psi   = 0;
+theta = 2.95;
+
+phi   = 0.1;
+
+P = [rho, tau; tau/2, rho/3];
+
+external_function(nargs=3, name=extFunNoDerivsMatrix, in_arg_dim=[2,2,2], out_arg_dim=2);
+
+model(use_dll);
+c*theta*h^(1+psi)=(1-alpha)*y;
+k = beta*(((exp(A[2])*c)/(exp(A[2](+1))*c(+1)))
+          *(exp(A[2](+1))*alpha*y(+1)+(1-delta)*k));
+y = exp(A[1])*(k(-1)^alpha)*(h^(1-alpha));
+k = exp(A[2])*(y-c)+(1-delta)*k(-1) ;
+A = extFunNoDerivsMatrix(P[:,1], transpose(P[:,2]), A(-1)) + E;
+end;
+
+initval;
+y = 1.08068253095672;
+c = 0.80359242014163;
+h = 0.29175631001732;
+k = 11.08360443260358;
+
+A[1] = 0;
+A[2] = 0;
+E[:] = 0;
+end;
+
+shocks;
+var E[:]; stderr 0.009;
+var E[1], E[2] = phi*0.009*0.009;
+end;
+
+stoch_simul;
+
+L = load('benchmark_results.mat');
+if max(max(abs(L.oo_.dr.ghu - oo_.dr.ghu))) > 1e-12
+  error('Failure in matrix notation with external function')
+end
+if max(max(abs(L.oo_.dr.ghx - oo_.dr.ghx))) > 1e-12
+  error('Failure in matrix notation with external function')
+end
+if max(max(abs(L.oo_.dr.ghxu - oo_.dr.ghxu))) > 1e-12
+  error('Failure in matrix notation with external function')
+end
+if max(max(abs(L.oo_.dr.ghxx - oo_.dr.ghxx))) > 1e-12
+  error('Failure in matrix notation with external function')
+end
+if max(max(abs(L.oo_.dr.ghuu - oo_.dr.ghuu))) > 1e-12
+  error('Failure in matrix notation with external function')
+end
+if max(max(abs(L.oo_.dr.ghs2 - oo_.dr.ghs2))) > 1e-12
+  error('Failure in matrix notation with external function')
+end
+if max(max(abs(L.oo_.var - oo_.var))) > 1e-12
+  error('Failure in matrix notation with external function')
+end
diff --git a/tests/matrix_notation/stochastic_matrix.mod b/tests/matrix_notation/stochastic_matrix.mod
new file mode 100644
index 0000000000000000000000000000000000000000..d055f5e60c4c5ee028b968792e63c76889589c61
--- /dev/null
+++ b/tests/matrix_notation/stochastic_matrix.mod
@@ -0,0 +1,46 @@
+// Do not forget to update stochastic_scalar.mod when this file is modified
+
+var(rows=3) X, Y, Z, T;
+var r;
+varexo(rows=3) E, U, V, W;
+varexo f;
+parameters(rows=3, cols=3) P;
+varexo_det(rows=3) D;
+
+P = [.1, .2, -.3; -.1, .5, .1; .2, -.3, .9];
+
+model;
+  X = P*X(-1)+E;
+  transpose(Y) = transpose(Y(-1))*transpose(P)+transpose(U);
+  Z[1] = V[1];
+  Z[2] = V[2];
+  Z[3] = V[3];
+  T = D + W;
+  r = 0.9*r(+1)+f;
+end;
+
+shocks;
+  var E = [ 0.1, 0.2, 0.3 ];
+  var E[1], E[2] = 0.08;
+
+  var U[1] = 0.4;
+  var U[2]; stderr 0.5;
+  var U[3]; stderr 0.6;
+  corr U[1], U[2] = 0.09;
+
+  var V; stderr [ 1.1, 1.2, 1.3 ];
+
+  vcov W = [
+      1, 0.2, 0.3;
+    0.2, 5,   0.6;
+    0.3, 0.6,   9
+  ];
+
+  var f = 3;
+  var f, E[2] = 0.5;
+  corr U[1], f = 0.7;
+end;
+
+stoch_simul(order=1, nodecomposition, irf=0);
+
+write_latex_dynamic_model;
diff --git a/tests/matrix_notation/stochastic_scalar.mod b/tests/matrix_notation/stochastic_scalar.mod
new file mode 100644
index 0000000000000000000000000000000000000000..3c5be1857c21938e5398bffb8f341eba5241d8bf
--- /dev/null
+++ b/tests/matrix_notation/stochastic_scalar.mod
@@ -0,0 +1,72 @@
+// Scalar translation of stochastic_matrix.mod
+// The declaration orders are the same as those of the expanded matrix version
+var r X_1 X_2 X_3 Y_1 Y_2 Y_3 Z_1 Z_2 Z_3 T_1 T_2 T_3;
+varexo f E_1 E_2 E_3 U_1 U_2 U_3 V_1 V_2 V_3 W_1 W_2 W_3;
+parameters P_1_1 P_1_2 P_1_3 P_2_1 P_2_2 P_2_3 P_3_1 P_3_2 P_3_3;
+varexo_det D_1 D_2 D_3;
+
+P_1_1 = .1;
+P_1_2 = .2;
+P_1_3 = -.3;
+P_2_1 = -.1;
+P_2_2 = .5;
+P_2_3 = .1;
+P_3_1 = .2;
+P_3_2 = -.3;
+P_3_3 = .9;
+
+model;
+  X_1 = P_1_1*X_1(-1) + P_1_2*X_2(-1) + P_1_3*X_3(-1) + E_1;
+  X_2 = P_2_1*X_1(-1) + P_2_2*X_2(-1) + P_2_3*X_3(-1) + E_2;
+  X_3 = P_3_1*X_1(-1) + P_3_2*X_2(-1) + P_3_3*X_3(-1) + E_3;
+  Y_1 = P_1_1*Y_1(-1) + P_1_2*Y_2(-1) + P_1_3*Y_3(-1) + U_1;
+  Y_2 = P_2_1*Y_1(-1) + P_2_2*Y_2(-1) + P_2_3*Y_3(-1) + U_2;
+  Y_3 = P_3_1*Y_1(-1) + P_3_2*Y_2(-1) + P_3_3*Y_3(-1) + U_3;
+  Z_1 = V_1;
+  Z_2 = V_2;
+  Z_3 = V_3;
+  T_1 = D_1 + W_1;
+  T_2 = D_2 + W_2;
+  T_3 = D_3 + W_3;
+  r = 0.9*r(+1)+f;
+end;
+
+shocks;
+  var E_1 = 0.1;
+  var E_2 = 0.2;
+  var E_3 = 0.3;
+  var E_1, E_2 = 0.08;
+
+  var U_1 = 0.4;
+  var U_2; stderr 0.5;
+  var U_3; stderr 0.6;
+  corr U_1, U_2 = 0.09;
+
+  var V_1; stderr 1.1;
+  var V_2; stderr 1.2;
+  var V_3; stderr 1.3;
+
+  var W_1 = 1;
+  var W_2 = 5;
+  var W_3 = 9;
+  var W_1, W_2 = 0.2;
+  var W_1, W_3 = 0.3;
+  var W_2, W_3 = 0.6;
+
+  var f = 3;
+  var f, E_2 = 0.5;
+  corr U_1, f = 0.7;
+end;
+
+stoch_simul(order=1, nodecomposition, irf=0);
+
+L = load('stochastic_matrix_results.mat');
+if max(max(abs(L.oo_.dr.ghu - oo_.dr.ghu))) > 1e-12
+  error('Failure in matrix expansion')
+end
+if max(max(abs(L.oo_.dr.ghx - oo_.dr.ghx))) > 1e-12
+  error('Failure in matrix expansion')
+end
+if max(max(abs(L.oo_.var - oo_.var))) > 1e-12
+  error('Failure in matrix expansion')
+end