diff --git a/test/test_analytical_derivatives.jl b/test/test_analytical_derivatives.jl
index 9fbf16ef..783e7a08 100644
--- a/test/test_analytical_derivatives.jl
+++ b/test/test_analytical_derivatives.jl
@@ -5,6 +5,7 @@ using Dynare.FastLapackInterface
 using Test
 using SparseArrays
 using SuiteSparse
+using GeneralizedSylvesterSolver
 
 #Model
 context = @dynare "models/analytical_derivatives/fs2000_sa.mod" "params_derivs_order=1" "notmpterms";
@@ -168,7 +169,6 @@ workspace2 = init_derivatives2_workspace(n, m)
 workspace3 = init_derivatives3_workspace(n, m)
 workspace4 = init_derivatives4_workspace(n, m)
 
-# Luego puedes utilizar workspace para llamar a la función Derivatives con tus matrices df_dx y df_dp
 @time sol1 = Derivatives(workspace1, df_dx, df_dp);
 @time sol2 = Derivatives2(workspace2, df_dx, df_dp);
 @time sol3 = Derivatives3(workspace3, df_dx, df_dp);
@@ -227,5 +227,55 @@ X = zeros(n, n)
 # set nonzero columns
 X[:, k1] .= context.results.model_results[1].linearrationalexpectations.g1_1
 
+#Generalized Sylvester: ax + bxc = d
+a = A*X + B
+b = Matrix(A)
+c = X
+X2 = X*X
+d = zeros(n, n, m)
+
+a_orig = copy(a)
+b_orig = copy(b)
+c_orig = copy(c)
+d_orig = copy(d)
+
+order=1
+ws = GeneralizedSylvesterWs(n,n,n,order)
+#Solve UQME using generalized_sylvester_solver!
+for i in 1:m
+    @views begin
+        mul!(d[:,:,i], dA_dp[:,:,i], X2)
+        mul!(d[:,:,i], dB_dp[:,:,i], X, true, true) 
+        d[:,:,i] .+= dC_dp[:,:,i] 
+    end
+    generalized_sylvester_solver!(a, b, c, d[:,:,i], order, ws)
+end
+
+#Test 1
+for i in 1:m
+    @test a_orig*d[:,:,i] + b_orig*d[:,:,i]*c_orig ≈ d_orig[:,:,i] #Fail
+end
+
+#Test 2
+for i in 1:m
+    @test d[:,:,i] ≈ reshape((kron(I(n^order),a_orig) + kron(c_orig',b_orig))\vec(d_orig[:,:,i]),n,n^order) #Fail
+end
+
+#Test 3
+using FiniteDifferences
+
+function funX(params)
+    X = zeros(n, n)
+    X[:, k1] .= context.results.model_results[1].linearrationalexpectations.g1_1
+    return X #vec()
+end
+
+fd = central_fdm(5, 1)
+dX_dz_tuple = jacobian(fd, funX, params)
+dX_dz_matrix = dX_dz_tuple[1]
+dX_dz = permutedims(reshape(dX_dz_matrix, m, n, n), (2, 3, 1))
+
+@test d ≈ dX_dz #Fail
+
 end # end module