diff --git a/code/m_chemistry.f90 b/code/m_chemistry.f90
index 4dd3dce..3d72f64 100644
--- a/code/m_chemistry.f90
+++ b/code/m_chemistry.f90
@@ -63,7 +63,8 @@ real function rate1_2step (ya, yx, theta)
     if(theta.lt.0.) t_reduce=0.
     if(theta.gt.1.) t_reduce=1.
 
-    rate1_2step = 0.
+    rate1_2step = lambda1 * y1 * y2 * & 
+                  exp (-(beta1*(1. - t_reduce))/(1. - hrp*(1. - t_reduce)))
 
 end function rate1_2step
 
@@ -84,7 +85,7 @@ real function rate2_2step (yx, theta)
     if(theta.lt.0.) t_reduce=0.
     if(theta.gt.1.) t_reduce=1.
 
-    rate2_2step = 0.
+    rate2_2step = lambda2 * yx * yx
 
 end function rate2_2step
 
diff --git a/code/m_parameters.f90 b/code/m_parameters.f90
index e8ce0d0..0367666 100644
--- a/code/m_parameters.f90
+++ b/code/m_parameters.f90
@@ -16,6 +16,8 @@ module m_parameters
 
      real :: pre,ac,bc,c_cut,c_ref
 
+     real :: lambda1, lambda2, beta1, hrp
+
      integer :: syp,eyp,twod
 
      real :: SL_u
@@ -28,7 +30,7 @@ module m_parameters
 
      REAL :: hx,dt,vis,sc,diff,tf,t_now,t_uf,dt_uf
 
-     REAL :: lewis
+     REAL :: lewis, le_a, le_x
 
      REAL :: minf,tar_lo,u0,ctmp,lo_flm=0.
 
@@ -44,6 +46,15 @@ module m_parameters
 
 contains
 
+     SUBROUTINE SET_CHEMISTRY
+
+       lambda1 = 100000.
+       lambda2 = lambda1 * 63.6
+       beta1 = 9.1
+       hrp = 0.8136765
+
+     END SUBROUTINE SET_CHEMISTRY
+
      SUBROUTINE READ_INTRO
        CHARACTER(LEN=8) :: cdum
        INTEGER :: itape=300,otape=301
@@ -96,6 +107,18 @@ contains
         WRITE(otape,*) cdum,prof_wr
         READ(itape,*) cdum,lewis
         WRITE(otape,*) cdum,lewis
+        READ(itape,*) cdum,lambda1
+        WRITE(otape,*) cdum,lambda1
+        READ(itape,*) cdum,lambda2
+        WRITE(otape,*) cdum,lambda2
+        READ(itape,*) cdum,beta1
+        WRITE(otape,*) cdum,beta1
+        READ(itape,*) cdum,hrp
+        WRITE(otape,*) cdum,hrp
+        READ(itape,*) cdum,le_a
+        WRITE(otape,*) cdum,le_a
+        READ(itape,*) cdum,le_x
+        WRITE(otape,*) cdum,le_x
 
        IF (.not.read_stdin) THEN
          CLOSE(itape)
diff --git a/code/sample/4pi-IC1/itape b/code/sample/4pi-IC1/itape
index 2c08701..7a31b89 100644
--- a/code/sample/4pi-IC1/itape
+++ b/code/sample/4pi-IC1/itape
@@ -9,7 +9,7 @@
         ac     26.7            !22.250          !22.25        !22.56        !22.3        !22.39
         bc     3.              !7.000           !7.0          !7.0          !7.0         !7.0
         u0     0.21
-        tf     80.
+        tf     200.
      dt_uf     0.4
       ctmp     0.        !1.0e-14
      c_cut     0.001 !0.012     ! c < c_cut -> wrate = 0.
@@ -17,4 +17,9 @@
     min_wr     0.        ! 5.0e-14
    prof_wr     1.0       
      lewis     1.0
-
+       lambda1   100000.
+       lambda2   63.6E+5
+       beta1   9.1
+       hrp   0.8136765
+       le_a  1.
+       le_x  0.2
diff --git a/code/ysolve.f90 b/code/ysolve.f90
index 3933233..01032e8 100644
--- a/code/ysolve.f90
+++ b/code/ysolve.f90
@@ -223,7 +223,7 @@
 
       OPEN (305,FILE='sfield.dat')
       DO i=1,nx
-       WRITE (305,'(e30.20)') y1(i,2)
+       WRITE (305,'(4e30.20)') i*hx, y1(i,1), y1(i,2), y1(i,3)
       ENDDO
       CLOSE (305)
 
@@ -273,14 +273,15 @@
      END SUBROUTINE write_pre
 
      SUBROUTINE SET_BC
-      y1(1,1)=1.-ctmp
-      y1(1,2)=1.
+      y1(1,1)=1.
+      y1(1,2)=0.
       y1(1,3)=0.
      END SUBROUTINE SET_BC
 
      SUBROUTINE SET_IC
       INTEGER :: i, ifl, si
       REAL :: xi
+      REAL :: max_ysum
 
       ! initial flame thickness (0.2 pi n grids)
       ifl=INT(2.0*pi*0.1/hx)
@@ -303,13 +304,20 @@
 
 !       y1(i,1,1)=(1.-xi)*1.
        y1(i,1)=(1.-xi)   ! reactant mass fraction
+
+       y1(i,2) = (1./2.) * (lambda1/lambda2) * y1(i,1) * &
+         exp (-(beta1*(1. - xi))/(1. - hrp*(1. - xi)))
+
+       y1(i,3) = xi
+
       ENDDO
 
-      ! initialize intermediate species field
-      y1(:,2) = y1(:,1)
-
-      ! initialize temperature field
-      y1(:,3) = 1.0 - y1(:,1)
+      ! initialize species field
+      ! max_ysum = maxval(y1(:,1) + y1(:,2))
+      ! y1(:,1) = y1(:,1) / maxval(y1(:,1) + y1(:,2))
+      ! y1(:,2) = y1(:,2) / maxval(y1(:,1) + y1(;,2))
+      ! y1(:,1) = y1(:,1) / (y1(1,1) + y1(1,2))
+      ! y1(:,2) = y1(:,2) / (y1(1,1) + y1(1,2))
 
       ! calculate flame position
       pflame=0.
@@ -423,21 +431,23 @@
 
       DO i=1,nx
         wrate=rate_1step(ux(1,i), ux(3,i))
+        wrate1=rate1_2step(ux(1,i), ux(2,i), ux(3,i))
+        wrate2=rate2_2step(ux(2,i), ux(3,i))
 
         ! - u*dY/dx + D*d2Y/d2x
-        f(i,1) = - ( u(i)*dux(1,i) ) + diff * d2ux(1,i) - wrate
+        f(i,1) = - ( u(i)*dux(1,i) ) + (diff/le_a) * d2ux(1,i) - wrate1
 
         ! - u*dY/dx + D*d2Y/d2x
-        f(i,2) = - ( u(i)*dux(2,i) ) + diff * d2ux(2,i) - wrate
+        f(i,2) = - ( u(i)*dux(2,i) ) + (diff/le_x) * d2ux(2,i) + wrate1 - 2.*wrate2
 
         ! -  u*dT/dx + D*d2T/d2x
-        f(i,3) = - ( u(i)*dux(3,i) ) + (lewis * diff) * d2ux(3,i) + wrate
+        f(i,3) = - ( u(i)*dux(3,i) ) + (diff) * d2ux(3,i) + 2.*wrate2
 
         ! 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)
-            f(nx,3) = wrate - u(nx)*dux(3,nx)
+            f(nx,1) = -wrate1 - u(nx)*dux(1,nx)
+            f(nx,2) = wrate1 - 2.*wrate2 - u(nx)*dux(2,nx)
+            f(nx,3) = 2.*wrate2 - u(nx)*dux(3,nx)
         ENDIF
       ENDDO