diff --git a/matlab/getPowerDeriv.m b/matlab/getPowerDeriv.m
index 3e0d1823e55a344c8e33236e4be8b93e0e681624..4dee1c4172694fea043a5b594f296bd8d130d846 100644
--- a/matlab/getPowerDeriv.m
+++ b/matlab/getPowerDeriv.m
@@ -30,7 +30,7 @@ function dxp=getPowerDeriv(x,p,k)
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
 
-if (abs(x) < 1e-12) && (p > 0) && (k > p) && (abs(p - round(p)) < 1e-12)
+if (abs(x) < 1e-12) && (p >= 0) && (k > p) && (abs(p - round(p)) < 1e-12)
     dxp = 0;
 else
     dxp = x^(p-k);
@@ -39,4 +39,23 @@ else
         p = p-1;
     end
 end
-end
+
+return % --*-- Unit tests --*--
+
+%@test:1
+x=getPowerDeriv(2,3,1);
+t(1)=all(abs(x-3*4)<1e-10)
+x=getPowerDeriv(0,2,2);
+t(2)=all(abs(x-2)<1e-10)
+x=getPowerDeriv(0,2,3); %special case evaluates to 0
+t(3)=all(abs(x-0)<1e-10)
+x=getPowerDeriv(1e-13,2,3-1e-13); %0 within tolerance
+t(4)=all(abs(x-0)<1e-10)
+x=getPowerDeriv(0,0,1);
+t(5)=all(abs(x-0)<1e-10)
+x=getPowerDeriv(0,0,0);
+t(6)=all(abs(x-1)<1e-10);
+x=getPowerDeriv(0,1/3,1); %derivative evaluating to Inf due to division by 0
+t(7)= isinf(x)
+T = all(t);
+%@eof:1
\ No newline at end of file