Index: trunk/src/wave_equations_integrator.F90
===================================================================
--- trunk/src/wave_equations_integrator.F90	(revision 585)
+++ trunk/src/wave_equations_integrator.F90	(revision 594)
@@ -345,4 +345,9 @@
     f(1) = r0**(s - M_ONE)*(b0 + b1*r0 + b2*r0**2)
 
+    ! MGGA term
+    if (wave_eq == SCHRODINGER) then
+      f(1) = f(1)*(M_ONE + M_TWO*vtau(potential, r0, qn))
+    end if
+
     !Set the parameters needed to compute the derivatives of the functions
     call set_derivs_params(qn, e, wave_eq, potential)
@@ -374,5 +379,7 @@
     case (SCHRODINGER)
       wf_out(1:nstep,1) = g(1:nstep)
-      wfp_out(1:nstep,1) = f(1:nstep)
+      do i = 1, nstep
+        wfp_out(i,1) = f(i)/(M_ONE + M_TWO*vtau(potential, r_out(i), qn))
+      end do
     case (SCALAR_REL)
       wf_out(1:nstep,1) = g(1:nstep)
@@ -439,5 +446,5 @@
     case (SCHRODINGER)
       g(1) = SCHRODINGER_FINF
-      f(1) = -sqrt(-M_TWO*e)*g(1)
+      f(1) = -sqrt(-M_TWO*e)*g(1)*(M_ONE + M_TWO*vtau(potential, rinf, qn))
     case (DIRAC)
       !Values of the wavefunctions at infinity
@@ -487,5 +494,7 @@
     case (SCHRODINGER)
       wf_in(1:nstep,1) = g(1:nstep)
-      wfp_in(1:nstep,1) = f(1:nstep)
+      do i = 1, nstep
+        wfp_in(i,1) = f(i)/(M_ONE + M_TWO*vtau(potential, r_in(i), qn))
+      end do
     case (SCALAR_REL)
       wf_in(1:nstep,1) = g(1:nstep)