diff --git a/src/polynomials.f90 b/src/polynomials.f90
index 1a023eba1..2426ba48d 100644
--- a/src/polynomials.f90
+++ b/src/polynomials.f90
@@ -141,7 +141,7 @@ subroutine selfTest()
   type(tPolynomial) :: p1, p2
   real(pReal), dimension(5) :: coef
   integer :: i
-  real(pReal) :: x_ref, x
+  real(pReal) :: x_ref, x, y
   class(tNode), pointer :: dict
   character(len=pStringLen), dimension(size(coef)) :: coef_s
   character(len=pStringLen) :: x_ref_s, x_s, YAML_s
@@ -155,6 +155,9 @@ subroutine selfTest()
   x_ref = x_ref*10_pReal -0.5_pReal
   x = x*10_pReal -0.5_pReal
 
+  p1 = polynomial([coef(1)],x_ref)
+  if (dNeq(p1%at(x),coef(1)))      error stop 'polynomial: eval(constant)'
+
   p1 = polynomial(coef,x_ref)
   if (dNeq(p1%at(x_ref),coef(1)))  error stop 'polynomial: @ref'
 
@@ -172,6 +175,8 @@ subroutine selfTest()
   Dict => YAML_parse_str(trim(YAML_s))
   p2 = polynomial(dict%asDict(),'C','T')
   if (dNeq(p1%at(x),p2%at(x),1.0e-10_pReal))                      error stop 'polynomials: init'
+  y = coef(1)+coef(2)*(x-x_ref)+coef(3)*(x-x_ref)**2+coef(4)*(x-x_ref)**3+coef(5)*(x-x_ref)**4
+  if (dNeq(p1%at(x),y,1.0e-7_pReal))                             error stop 'polynomials: eval(full)'
 
   YAML_s = 'C: 0.0'//IO_EOL//&
            'C,T: '//trim(adjustl(coef_s(2)))//IO_EOL//&
@@ -189,15 +194,15 @@ subroutine selfTest()
   if (dNeq(p1%at(x_ref+x),p1%at(x_ref-x),1e-10_pReal))            error stop 'polynomials: eval(quadratic)'
   if (dNeq(p1%der1_at(x_ref+x),-p1%der1_at(x_ref-x),1e-10_pReal)) error stop 'polynomials: eval_der(quadratic)'
 
-  YAML_s = 'Y: 0.0'//IO_EOL//&
-           'Y,X^3: '//trim(adjustl(coef_s(1)))//IO_EOL//&
+  YAML_s = 'Y: '//trim(adjustl(coef_s(1)))//IO_EOL//&
+           'Y,X^3: '//trim(adjustl(coef_s(2)))//IO_EOL//&
            'X_ref: '//trim(adjustl(x_ref_s))//IO_EOL
   Dict => YAML_parse_str(trim(YAML_s))
   p1 = polynomial(dict%asDict(),'Y','X')
-  if (dNeq(p1%at(x_ref+x),-p1%at(x_ref-x),1.0e-8_pReal))          error stop 'polynomials: eval(cubic)'
+  if (dNeq(p1%at(x_ref+x)-coef(1),-(p1%at(x_ref-x)-coef(1)),1.0e-8_pReal))          error stop 'polynomials: eval(cubic)'
 
-  YAML_s = 'Y: 0.0'//IO_EOL//&
-           'Y,X^4: '//trim(adjustl(coef_s(1)))//IO_EOL//&
+  YAML_s = 'Y: '//trim(adjustl(coef_s(1)))//IO_EOL//&
+           'Y,X^4: '//trim(adjustl(coef_s(2)))//IO_EOL//&
            'X_ref: '//trim(adjustl(x_ref_s))//IO_EOL
   Dict => YAML_parse_str(trim(YAML_s))
   p1 = polynomial(dict%asDict(),'Y','X')