diff --git a/code/ysolve.f90 b/code/ysolve.f90
index 5aee200..25ea98f 100644
--- a/code/ysolve.f90
+++ b/code/ysolve.f90
@@ -274,6 +274,7 @@
 
      SUBROUTINE SET_BC
       y1(1,1)=1.-ctmp
+      y1(1,2)=0.
      END SUBROUTINE SET_BC
 
      SUBROUTINE SET_IC
@@ -406,45 +407,36 @@
       REAL :: wrate1, wrate2
 
 ! x-direction
-        DO i=1,nx
+      DO i=1,nx
          ux(1,i)=r1(i,1)                              ! Y
          ux(2,i)=r1(i,2)                              ! T
          ux(3,i)=r1(i,3)                              ! T
-        ENDDO
+      ENDDO
 
-        CALL dfnonp(nx,hx,ux(:,:),dux(:,:),3,1)
+      CALL dfnonp(nx,hx,ux(:,:),dux(:,:),3,1)
 
-        CALL d2fnonp(nx,hx,ux(:,:),d2ux(:,:),3,1)
+      CALL d2fnonp(nx,hx,ux(:,:),d2ux(:,:),3,1)
 
-        DO i=1,nx
-         wrate=rate_1step(ux(1,i), ux(2,i))
+      DO i=1,nx
+        wrate=rate_1step(ux(1,i), ux(2,i))
 
-! -0.5*( d(u*Y)/dx + u*dY/dx + Y*du/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,1) = - ( u(i)*dux(1,i) ) + diff * d2ux(1,i) - wrate
 
-! -  u*dT/dx + D*d2T/d2x
-         f(i,2) = - ( u(i)*dux(2,i) ) + (lewis * diff) * d2ux(2,i) + wrate
+        ! -  u*dT/dx + D*d2T/d2x
+        f(i,2) = - ( u(i)*dux(2,i) ) + (lewis * diff) * d2ux(2,i) + wrate
 
-         IF (i.eq.nx) THEN
-          Ly=u(nx)*dux(1,nx)                           ! Ly = u*dYr/dx
-          f(nx,1)=-wrate
+        ! Boundary conditions
+        IF (i.eq.nx) THEN
+            f(nx,1) = -wrate - u(nx)*dux(1,nx)
+            f(nx,2) = wrate - u(nx)*dux(2,nx)
+        ENDIF
+      ENDDO
 
-          Lt=u(nx)*dux(2,nx)                           ! Ly = u*dYr/dx
-          f(nx,2)=wrate
-         ENDIF
-        ENDDO
-
-!    Boundary conditionS
+      ! Boundary conditions
       f(1,1)=0.
-
-      Dy=Ly
-      f(nx,1)=f(nx,1)-Dy
-
-
       f(1,2)=0.
-
-      Dt=Lt
-      f(nx,2)=f(nx,2)-Dt
+      f(1,3)=0.
 
      END SUBROUTINE fns