diff --git a/code/ysolve.f90 b/code/ysolve.f90
index f066bcd..be1c76d 100644
--- a/code/ysolve.f90
+++ b/code/ysolve.f90
@@ -463,6 +463,46 @@
 
      END SUBROUTINE RK4
 !------------------------------------------------------------------------
+     SUBROUTINE fonestep(r1,f)
+      REAL, INTENT(IN),DIMENSION(:,:) :: r1
+      REAL, INTENT(OUT),DIMENSION(:,:) :: f
+      REAL, DIMENSION(nsp,nx) :: ux, dux, d2ux
+      INTEGER :: i,j,k
+      REAL :: wrate,Ly,Dy,Lt,Dt
+      REAL :: wrate1, wrate2
+
+! x-direction
+      DO i=1,nx
+         ux(1,i)=r1(i,1)                              ! Y
+         ux(2,i)=r1(i,2)                              ! T
+      ENDDO
+
+      CALL dfnonp(nx,hx,ux(:,:),dux(:,:),nsp,1)
+
+      CALL d2fnonp(nx,hx,ux(:,:),d2ux(:,:),nsp,1)
+
+      DO i=1,nx
+        wrate=rate_1step(ux(1,i), ux(2,i))
+
+        ! - u*dY/dx + D*d2Y/d2x
+        f(i,1) = - ( u(i)*dux(1,i) ) + (diff) * d2ux(1,i) - wrate
+
+        ! - u*dY/dx + D*d2Y/d2x
+        f(i,2) = - ( u(i)*dux(2,i) ) + (diff) * d2ux(2,i) + wrate
+
+        ! Boundary conditions
+        IF (i.eq.nx) THEN
+            f(nx,1) = -wrate1 - u(nx)*dux(1,nx)
+            f(nx,2) = wrate1 - u(nx)*dux(2,nx)
+        ENDIF
+      ENDDO
+
+      ! Boundary conditions
+      f(1,1)=0.
+      f(1,2)=0.
+
+     END SUBROUTINE fonestep
+
      SUBROUTINE fns(r1,f)
       REAL, INTENT(IN),DIMENSION(:,:) :: r1
       REAL, INTENT(OUT),DIMENSION(:,:) :: f