diff --git a/reactors/1D_pfr_surfchem.ipynb b/reactors/1D_pfr_surfchem.ipynb
index 75a8b57..09b1e9a 100644
--- a/reactors/1D_pfr_surfchem.ipynb
+++ b/reactors/1D_pfr_surfchem.ipynb
@@ -36,8 +36,6 @@
     "import cantera as ct\n",
     "import matplotlib.pyplot as plt\n",
     "%matplotlib inline\n",
-    "from prettyplotlib import brewer2mpl\n",
-    "colors = brewer2mpl.get_map('Set2','qualitative',8).mpl_colors\n",
     "print('Runnning Cantera version: ' + ct.__version__)"
    ]
   },
@@ -364,7 +362,8 @@
     "    old_api=False  # Forces use of new api (namedtuple)\n",
     ")\n",
     "\n",
-    "solution = solver.solve(np.arange(0,0.7,0.1), vec0,vecp0)"
+    "times = np.arange(0,0.7,0.01)\n",
+    "solution = solver.solve(times, vec0, vecp0)"
    ]
   },
   {
@@ -393,23 +392,17 @@
     }
    ],
    "source": [
-    "plt.rcParams['axes.labelsize'] = 18\n",
-    "plt.rcParams['xtick.labelsize'] = 16\n",
-    "plt.rcParams['ytick.labelsize'] = 16\n",
-    "\n",
     "# plot velocity of gas along the flow direction\n",
-    "plt.figure(figsize = (16,18))\n",
-    "plt.subplot(3,2,1)\n",
-    "plt.plot(np.arange(0,0.7,0.1), solution.values.y[:,0], color=colors[0])\n",
-    "plt.xlabel('Distance (m)')\n",
-    "plt.ylabel('Velocity (m/s)')\n",
+    "f, ax = plt.subplots(3,2, figsize=(9,9), dpi=96)\n",
+    "ax[0,0].plot(times, solution.values.y[:,0], color='C0')\n",
+    "ax[0,0].set_xlabel('Distance (m)')\n",
+    "ax[0,0].set_ylabel('Velocity (m/s)')\n",
     "\n",
     "# plot gas density along the flow direction\n",
-    "plt.subplot(3,2,2)\n",
-    "plt.plot(np.arange(0,0.7,0.1),solution.values.y[:,1], color=colors[1])\n",
-    "plt.xlabel('Distance (m)')\n",
-    "plt.ylabel('Density ($\\mathregular{kg/m^3}$)')\n",
-    "plt.ticklabel_format(axis='y', style='sci', scilimits=(-2,2))  # scientific notation\n",
+    "ax[0,1].plot(times, solution.values.y[:,1], color='C1')\n",
+    "ax[0,1].set_xlabel('Distance (m)')\n",
+    "ax[0,1].set_ylabel('Density ($\\mathregular{kg/m^3}$)')\n",
+    "ax[0,1].ticklabel_format(axis='y', style='sci', scilimits=(-2,2))  # scientific notation\n",
     "\n",
     "# plot major and minor gas species separately\n",
     "minor_idx = []\n",
@@ -422,36 +415,33 @@
     "        major_idx.append(i)\n",
     "\n",
     "# plot minor species\n",
-    "plt.subplot(3,2,3)\n",
     "for i in minor_idx:\n",
-    "    plt.plot(np.arange(0,0.7,0.1), solution.values.y[:,2+i], label=gas.species_names[i])\n",
-    "plt.legend(fontsize=11, loc='upper left')\n",
-    "plt.xlabel('Distance (m)')\n",
-    "plt.ylabel('Mass Fraction')\n",
-    "plt.ticklabel_format(axis='y', style='sci', scilimits=(-2,2))  # scientific notation\n",
+    "    style = '-' if i < 10 else '--'    \n",
+    "    ax[1,0].plot(times, solution.values.y[:,2+i], label=gas.species_names[i], linestyle=style)\n",
+    "ax[1,0].legend(fontsize=7, loc='upper right')\n",
+    "ax[1,0].set_xlabel('Distance (m)')\n",
+    "ax[1,0].set_ylabel('Mass Fraction')\n",
+    "ax[1,0].ticklabel_format(axis='y', style='sci', scilimits=(-2,2))  # scientific notation\n",
     "\n",
     "# plot major species\n",
-    "plt.subplot(3,2,4)\n",
     "for j in major_idx:\n",
-    "    plt.plot(np.arange(0,0.7,0.1),solution.values.y[:,2+j], label=gas.species_names[j])\n",
-    "plt.legend(fontsize=11, loc='best')\n",
-    "plt.xlabel('Distance (m)')\n",
-    "plt.ylabel('Mass Fraction')\n",
+    "    ax[1,1].plot(times,solution.values.y[:,2+j], label=gas.species_names[j])\n",
+    "ax[1,1].legend(loc='best')\n",
+    "ax[1,1].set_xlabel('Distance (m)')\n",
+    "ax[1,1].set_ylabel('Mass Fraction')\n",
     "\n",
     "# plot the pressure of the gas along the flow direction\n",
-    "plt.subplot(3,2,5)\n",
-    "plt.plot(np.arange(0,0.7,0.1), solution.values.y[:,2+N], color=colors[2])\n",
-    "plt.xlabel('Distance (m)')\n",
-    "plt.ylabel('Pressure (Pa)')\n",
+    "ax[2,0].plot(times, solution.values.y[:,2+N], color='C2')\n",
+    "ax[2,0].set_xlabel('Distance (m)')\n",
+    "ax[2,0].set_ylabel('Pressure (Pa)')\n",
     "\n",
     "# plot the site fraction of the surface species along the flow direction \n",
-    "plt.subplot(3,2,6)\n",
     "for i,name in enumerate(gas_Si_N_interface.species_names):\n",
-    "    plt.plot(np.arange(0,0.7,0.1), solution.values.y[:,3+N+i], label= '%s'%(name))\n",
-    "plt.legend(fontsize=12)\n",
-    "plt.xlabel('Distance (m)')\n",
-    "plt.ylabel('Site Fraction')\n",
-    "plt.show()"
+    "    ax[2,1].plot(times, solution.values.y[:,3+N+i], label=name)\n",
+    "ax[2,1].legend()\n",
+    "ax[2,1].set_xlabel('Distance (m)')\n",
+    "ax[2,1].set_ylabel('Site Fraction')\n",
+    "f.tight_layout(pad=0.5)"
    ]
   },
   {
@@ -672,15 +662,24 @@
     "solver = dae(\n",
     "    'ida',\n",
     "    residual, \n",
-    "    first_step_size=1e-18,\n",
     "    atol=1e-8,  # absolute tolerance for solution\n",
     "    rtol=1e-8,  # relative tolerance for solution\n",
     "    algebraic_vars_idx=[np.arange(3+N,3+N+M,1)], \n",
     "    max_steps=5000,\n",
+    "    one_step_compute=True,\n",
     "    old_api=False\n",
     ")\n",
     "\n",
-    "solution = solver.solve(np.arange(0,0.7,0.1), vec0,vecp0)"
+    "time = []\n",
+    "solution = []\n",
+    "state = solver.init_step(0.0, vec0, vecp0)\n",
+    "while state.values.t < 0.7:\n",
+    "    time.append(state.values.t)\n",
+    "    solution.append(state.values.y)\n",
+    "    state = solver.step(0.7)\n",
+    "\n",
+    "time = np.array(time)\n",
+    "solution = np.array(solution)"
    ]
   },
   {
@@ -703,71 +702,62 @@
     }
    ],
    "source": [
-    "plt.rcParams['axes.labelsize'] = 18\n",
-    "plt.rcParams['xtick.labelsize'] = 16\n",
-    "plt.rcParams['ytick.labelsize'] = 16\n",
+    "f, ax = plt.subplots(4,2, figsize=(9,12), dpi=96)\n",
     "\n",
     "# plot gas velocity along the flow direction\n",
-    "plt.figure(figsize = (16,24))\n",
-    "plt.subplot(4,2,1)\n",
-    "plt.plot(np.arange(0,0.7,0.1),solution.values.y[:,0], color=colors[0])\n",
-    "plt.xlabel('Distance (m)')\n",
-    "plt.ylabel('Velocity (m/s)')\n",
+    "ax[0,0].plot(time, solution[:,0], color='C0')\n",
+    "ax[0,0].set_xlabel('Distance (m)')\n",
+    "ax[0,0].set_ylabel('Velocity (m/s)')\n",
     "\n",
     "# plot gas density along the flow direction\n",
-    "plt.subplot(4,2,2)\n",
-    "plt.plot(np.arange(0,0.7,0.1),solution.values.y[:,1], color=colors[1])\n",
-    "plt.xlabel('Distance (m)')\n",
-    "plt.ylabel('Density ($\\mathregular{kg/m^3}$)')\n",
-    "plt.ticklabel_format(axis='y', style='sci', scilimits=(-2,2))  # scientific notation\n",
+    "ax[0,1].plot(time, solution[:,1], color='C1')\n",
+    "ax[0,1].set_xlabel('Distance (m)')\n",
+    "ax[0,1].set_ylabel('Density ($\\mathregular{kg/m^3}$)')\n",
+    "ax[0,1].ticklabel_format(axis='y', style='sci', scilimits=(-2,2))  # scientific notation\n",
     "\n",
     "# plot major and minor gas species separately\n",
     "minor_idx = []\n",
     "major_idx = []\n",
     "for i,name in enumerate(gas.species_names): \n",
-    "    mean = np.mean(solution.values.y[:,2+i])\n",
+    "    mean = np.mean(solution[:,2+i])\n",
     "    if mean <= 0.01:\n",
     "        minor_idx.append(i) \n",
     "    else:\n",
     "        major_idx.append(i)\n",
     "\n",
     "# plot minor gas species along the flow direction\n",
-    "plt.subplot(4,2,3)\n",
     "for i in minor_idx:\n",
-    "    plt.plot(np.arange(0,0.7,0.1), solution.values.y[:,2+i], label=gas.species_names[i])\n",
-    "plt.legend(fontsize=11, loc='best')\n",
-    "plt.xlabel('Distance (m)')\n",
-    "plt.ylabel('Mass Fraction')\n",
-    "plt.ticklabel_format(axis='y', style='sci', scilimits=(-2,2))  # scientific notation\n",
+    "    style = '-' if i < 10 else '--'\n",
+    "    ax[1,0].plot(time, solution[:,2+i], label=gas.species_names[i], linestyle=style)\n",
+    "ax[1,0].legend(fontsize=7.5, loc='best')\n",
+    "ax[1,0].set_xlabel('Distance (m)')\n",
+    "ax[1,0].set_ylabel('Mass Fraction')\n",
+    "ax[1,0].ticklabel_format(axis='y', style='sci', scilimits=(-2,2))  # scientific notation\n",
     "\n",
     "# plot major gas species along the flow direction\n",
-    "plt.subplot(4,2,4)\n",
     "for j in major_idx:\n",
-    "    plt.plot(np.arange(0,0.7,0.1), solution.values.y[:,2+j], label=gas.species_names[j])\n",
-    "plt.legend(fontsize=11, loc='best')\n",
-    "plt.xlabel('Distance (m)')\n",
-    "plt.ylabel('Mass Fraction')\n",
+    "    ax[1,1].plot(time, solution[:,2+j], label=gas.species_names[j])\n",
+    "ax[1,1].legend(fontsize=8, loc='best')\n",
+    "ax[1,1].set_xlabel('Distance (m)')\n",
+    "ax[1,1].set_ylabel('Mass Fraction')\n",
     "\n",
     "# plot the pressure of the gas along the flow direction\n",
-    "plt.subplot(4,2,5)\n",
-    "plt.plot(np.arange(0,0.7,0.1), solution.values.y[:,2+N], color=colors[2])\n",
-    "plt.xlabel('Distance (m)')\n",
-    "plt.ylabel('Pressure (Pa)')\n",
+    "ax[2,0].plot(time, solution[:,2+N], color='C2')\n",
+    "ax[2,0].set_xlabel('Distance (m)')\n",
+    "ax[2,0].set_ylabel('Pressure (Pa)')\n",
     "\n",
     "# plot the site fraction of the surface species along the flow direction\n",
-    "plt.subplot(4,2,6)\n",
     "for i,name in enumerate(gas_Si_N_interface.species_names):\n",
-    "    plt.plot(np.arange(0,0.7,0.1), solution.values.y[:,3+N+i], label=name)\n",
-    "plt.legend(fontsize=12)\n",
-    "plt.xlabel('Distance (m)')\n",
-    "plt.ylabel('Site Fraction')\n",
+    "    ax[2,1].plot(time, solution[:,3+N+i], label=name)\n",
+    "ax[2,1].legend(fontsize=8)\n",
+    "ax[2,1].set_xlabel('Distance (m)')\n",
+    "ax[2,1].set_ylabel('Site Fraction')\n",
     "\n",
     "# plot the temperature profile along the flow direction\n",
-    "plt.subplot(4,2,7)\n",
-    "plt.plot(np.arange(0,0.7,0.1), solution.values.y[:,-1], color=colors[3])\n",
-    "plt.xlabel('Distance (m)')\n",
-    "plt.ylabel('Temp (K)')\n",
-    "plt.show()"
+    "ax[3,0].plot(time, solution[:,-1], color='C3')\n",
+    "ax[3,0].set_xlabel('Distance (m)')\n",
+    "ax[3,0].set_ylabel('Temperature (K)')\n",
+    "f.tight_layout(pad=0.5)"
    ]
   }
  ],