Figure 5.4 (page 199):

Reaction-coordinate diagram.

Code for Figure 5.4

Text of the GNU GPL.

main.m

%% Copyright (C) 2001, James B. Rawlings and John G. Ekerdt
%%
%% This program is free software; you can redistribute it and/or
%% modify it under the terms of the GNU General Public License as
%% published by the Free Software Foundation; either version 2, or (at
%% your option) any later version.
%%
%% This program is distributed in the hope that it will be useful, but
%% WITHOUT ANY WARRANTY; without even the implied warranty of
%% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
%% General Public License for more details.
%%
%% You should have received a copy of the GNU General Public License
%% along with this program; see the file COPYING.  If not, write to
%% the Free Software Foundation, 59 Temple Place - Suite 330, Boston,
%% MA 02111-1307, USA.

initdata ();

en1 = 3;
en2 = 5;

[s, V, d3trace, num] = inittrace (en1, en2);

%----------------------------
				%the following part is different from
				%morse.m, hfc.m and hfs.m
hfptable = [s, V];
				%End of different part
%----------------------------

initdata.m

%% Copyright (C) 2001, James B. Rawlings and John G. Ekerdt
%%
%% This program is free software; you can redistribute it and/or
%% modify it under the terms of the GNU General Public License as
%% published by the Free Software Foundation; either version 2, or (at
%% your option) any later version.
%%
%% This program is distributed in the hope that it will be useful, but
%% WITHOUT ANY WARRANTY; without even the implied warranty of
%% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
%% General Public License for more details.
%%
%% You should have received a copy of the GNU General Public License
%% along with this program; see the file COPYING.  If not, write to
%% the Free Software Foundation, 59 Temple Place - Suite 330, Boston,
%% MA 02111-1307, USA.

function initdata ()

  global D1 D3 B1 B3 r10 r30 R_1 C Si A

  D1 = [129.3; 94.8; 129.3];
  B1 = [2.316; 2.047; 2.316];
  r10 = [0.931; 0.753; 0.931];
  D3 = [51.83; 34.56; 51.83];
  B3 = [2.229; 1.522; 2.229];
  r30 = [0.931; 0.753; 0.931];
  R_1 = [1.217; 0.443; 1.217];
  C = [1174.9; 388.8; 1174.9];
  Si = [3.090; 1.929; 3.090];
  A = [1.315; 0.762; 1.315];

inittrace.m

%% Copyright (C) 2001, James B. Rawlings and John G. Ekerdt
%%
%% This program is free software; you can redistribute it and/or
%% modify it under the terms of the GNU General Public License as
%% published by the Free Software Foundation; either version 2, or (at
%% your option) any later version.
%%
%% This program is distributed in the hope that it will be useful, but
%% WITHOUT ANY WARRANTY; without even the implied warranty of
%% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
%% General Public License for more details.
%%
%% You should have received a copy of the GNU General Public License
%% along with this program; see the file COPYING.  If not, write to
%% the Free Software Foundation, 59 Temple Place - Suite 330, Boston,
%% MA 02111-1307, USA.

function [s, V, d3trace, num] = inittrace (en1, en2)

  global D1 D3 B1 B3 r10 r30 R_1 C Si A

  R1sar = 1.29781;
  R2sar = 0.80363;
  numr = 30;

  R1r = linspace (R1sar, en1, numr)';

  opts = odeset ('AbsTol', sqrt (eps), 'RelTol', sqrt (eps));
  [tsolver, R2r] = ode15s (@rder, R1r, R2sar, opts);

  R1sal = 1.297;
  R2sal = 0.80363;
  numl = 1000;

  R2l = linspace (R2sal, en2, numl)';

  opts  =  odeset ('AbsTol', sqrt (eps), 'RelTol', sqrt (eps));
  [tsolver, R1l] = ode15s (@lder, R2l, R1sal, opts);

  R1rn = zeros (numr, 1);
  R2rn = zeros (numr, 1);
  num = numl+numr;

  for i = 1:numr
    R1rn(i) = R1r(numr+1-i);
    R2rn(i) = R2r(numr+1-i);
  end

  R1 = [R1rn; R1l];
  R2 = [R2rn; R2l];
  trtable = [R1, R2];

  s = zeros (num, 1);
  V = zeros (num, 1);
  D = zeros (num, 1);

  s(1) = 0;
  V(1) = Vau (R1(1), R2(1));
  D(1) = 1;

  for i = 2:num
    D(i) = i;
    s(i) = sqrt ((R1(i)-R1(i-1))^2 + (R2(i)-R2(i-1))^2) + s(i-1);
    V(i) = Vau (R1(i), R2(i));
  end

  pot = [s, V, R1, R2, D];

  for i = 1:num
    j = i*3-2;
    d3trace(j) = R1(i);
    d3trace(j+1) = R2(i);
    d3trace(j+2) = V(i);
  end

  d3trace = d3trace';

  ## this seems to be unused.
  ## for a = 1:3
  ##   for i = 1:num
  ##    j = i*3-2;
  ##    d3sptrace(a,j) = R1(i)+(a-2)*0.1;
  ##    d3sptrace(a,j+1) = R2(i);
  ##    d3sptrace(a,j+2) = -100+(a-2);
  ##   end
  ## end

lder.m

%% Copyright (C) 2001, James B. Rawlings and John G. Ekerdt
%%
%% This program is free software; you can redistribute it and/or
%% modify it under the terms of the GNU General Public License as
%% published by the Free Software Foundation; either version 2, or (at
%% your option) any later version.
%%
%% This program is distributed in the hope that it will be useful, but
%% WITHOUT ANY WARRANTY; without even the implied warranty of
%% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
%% General Public License for more details.
%%
%% You should have received a copy of the GNU General Public License
%% along with this program; see the file COPYING.  If not, write to
%% the Free Software Foundation, 59 Temple Place - Suite 330, Boston,
%% MA 02111-1307, USA.

function retval = lder (r2, r1)

  global D1 D3 B1 B3 r10 r30 R_1 C Si A

  r3 = r1+r2;
  r(1) = r1;
  r(2) = r2;
  r(3) = r3;

  for x = 1:3
   E1(x) = D1(x)*(exp(-2*B1(x)*(r(x)-r10(x)))-2*exp(-B1(x)*(r(x)-r10(x))));
   E3(x,1) = D3(x)*(exp(-2*B3(x)*(r(x)-r30(x)))+2*exp(-B3(x)*(r(x)-r30(x))));
   E3(x,2) = C(x)*(r(x)+A(x))*exp(-Si(x)*r(x));
    if(r(x) <= R_1(x))
     y = 1;
     else 
     y = 2;
    end
   Q(x) = (E1(x)+E3(x,y))/2;
   J(x) = (E1(x)-E3(x,y))/2;
  end

  %dQ(a,b) means dQ(a)/dr(b)
  for i = 1:3
     dE1(i,i) = -2*B1(i)*D1(i)*(exp(-2*B1(i)*(r(i)-r10(i)))-exp(-B1(i)*(r(i)-r10(i))));
   if(r(i) <= R_1(i))
     dE3(i,i) = -2*B3(i)*D3(i)*(exp(-2*B3(i)*(r(i)-r30(i)))+exp(-B3(i)*(r(i)-r30(i))));
   else
     dE3(i,i) = C(i)*exp(-Si(i)*r(i))*(1-Si(i)*(r(i)+A(i)));
   end
  end
  dE1(1,2) = 0;
  dE3(1,2) = 0;
  dE1(2,1) = 0;
  dE3(2,1) = 0;
  dE1(3,1) = dE1(3,3);
  dE3(3,1) = dE3(3,3);
  dE1(3,2) = dE1(3,3);
  dE3(3,2) = dE3(3,3);
  dE1(1,3) = 1000;
  dE3(1,3) = 1000;
  dE1(2,3) = 1000;
  dE3(2,3) = 1000;
  dE1(3,3) = 1000;
  dE3(3,3) = 1000;
  for i = 1:3
   for j = 1:3
    dQ(i,j) = (dE1(i,j)+dE3(i,j))/2;
    dJ(i,j) = (dE1(i,j)-dE3(i,j))/2;
   end
  end

  dV(1) = dQ(1,1)+dQ(3,1)-0.5/sqrt(0.5*((J(1)-J(2))*(J(1)-J(2))+(J(2)-J(3))*(J(2)-J(3))+(J(1)-J(3))*(J(1)-J(3))))*(J(1)*(2*dJ(1,1)-dJ(3,1))+J(2)*(-dJ(1,1)-dJ(3,1))+J(3)*(2*dJ(3,1)-dJ(1,1)));

  dV(2) = dQ(2,2)+dQ(3,2)-0.5/sqrt(0.5*((J(1)-J(2))*(J(1)-J(2))+(J(2)-J(3))*(J(2)-J(3))+(J(1)-J(3))*(J(1)-J(3))))*(J(1)*(-dJ(2,2)-dJ(3,2))+J(2)*(2*dJ(2,2)-dJ(3,2))+J(3)*(2*dJ(3,2)-dJ(2,2)));

  retval = dV(1)/dV(2);

rder.m

%% Copyright (C) 2001, James B. Rawlings and John G. Ekerdt
%%
%% This program is free software; you can redistribute it and/or
%% modify it under the terms of the GNU General Public License as
%% published by the Free Software Foundation; either version 2, or (at
%% your option) any later version.
%%
%% This program is distributed in the hope that it will be useful, but
%% WITHOUT ANY WARRANTY; without even the implied warranty of
%% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
%% General Public License for more details.
%%
%% You should have received a copy of the GNU General Public License
%% along with this program; see the file COPYING.  If not, write to
%% the Free Software Foundation, 59 Temple Place - Suite 330, Boston,
%% MA 02111-1307, USA.

function retval = rder (r1, r2)

  global D1 D3 B1 B3 r10 r30 R_1 C Si A

  retval = 1 / lder (r2, r1, D1, D3, B1, B3, r10, r30, R_1, C, Si, A);

Matlab dynamically-linked MEX functions (C language)

Vau.c

/*

Copyright (C) 2001, James B. Rawlings and John G. Ekerdt

This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License as
published by the Free Software Foundation; either version 2, or (at
your option) any later version.

This program is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program; see the file COPYING.  If not, write to
the Free Software Foundation, 59 Temple Place - Suite 330, Boston,
MA 02111-1307, USA.

*/

#include 

#include 

void
mexFunction (int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[])
{
  int nargin = nrhs;

  if (nargin == 2)
    {
      double r1 = mxGetScalar (prhs[0]);
      double r2 = mxGetScalar (prhs[1]);

      mxArray *mx_D1 = mexGetVariable ("global", "D1");
      mxArray *mx_D3 = mexGetVariable ("global", "D3");
      mxArray *mx_B1 = mexGetVariable ("global", "B1");
      mxArray *mx_B3 = mexGetVariable ("global", "B3");
      mxArray *mx_r10 = mexGetVariable ("global", "r10");
      mxArray *mx_r30 = mexGetVariable ("global", "r30");
      mxArray *mx_R_1 = mexGetVariable ("global", "R_1");
      mxArray *mx_C = mexGetVariable ("global", "C");
      mxArray *mx_Si = mexGetVariable ("global", "Si");
      mxArray *mx_A = mexGetVariable ("global", "A");

      double *D1 = mxGetPr (mx_D1);
      double *D3 = mxGetPr (mx_D3);
      double *B1 = mxGetPr (mx_B1);
      double *B3 = mxGetPr (mx_B3);
      double *r10 = mxGetPr (mx_r10);
      double *r30 = mxGetPr (mx_r30);
      double *R_1 = mxGetPr (mx_R_1);
      double *C = mxGetPr (mx_C);
      double *Si = mxGetPr (mx_Si);
      double *A = mxGetPr (mx_A);

      double r3 = r1 + r2;

      double r[3];

      double E1[3];
      double E3[3][2];
      double Q[3];
      double J[3];

      int i;

      double retval;

      r[0] = r1;
      r[1] = r2;
      r[2] = r3;

      for (i = 0; i < 3; i++)
	{
	  int j;

	  E1[i] = D1[i]*(exp(-2*B1[i]*(r[i]-r10[i]))-2*exp(-B1[i]*(r[i]-r10[i])));
	  E3[i][0] = D3[i]*(exp(-2*B3[i]*(r[i]-r30[i]))+2*exp(-B3[i]*(r[i]-r30[i])));
	  E3[i][1] = C[i]*(r[i]+A[i])*exp(-Si[i]*r[i]);

	  j = (r[i] > R_1[i]);

	  Q[i] = (E1[i] + E3[i][j]) / 2;
	  J[i] = (E1[i] - E3[i][j]) / 2;
	}

      retval = Q[0]+Q[1]+Q[2]-sqrt(0.5*((J[0]-J[1])*(J[0]-J[1])+(J[1]-J[2])*(J[1]-J[2])+(J[0]-J[2])*(J[0]-J[2])));

      plhs[0] = mxCreateDoubleScalar (retval);
    }
  else
    mexErrMsgTxt ("usage: Vau (r1, r2)");
}

Octave dynamically-linked functions (C++ language)

Vau.cc

/*

Copyright (C) 2001, James B. Rawlings and John G. Ekerdt

This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License as
published by the Free Software Foundation; either version 2, or (at
your option) any later version.

This program is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program; see the file COPYING.  If not, write to
the Free Software Foundation, 59 Temple Place - Suite 330, Boston,
MA 02111-1307, USA.

*/

#include 
#include 

DEFUN_DLD (Vau, args, ,
  "")
{
  octave_value retval;

  int nargin = args.length ();

  if (nargin == 2)
    {
      double r1 = args(0).double_value ();
      double r2 = args(1).double_value ();

      if (! error_state)
	{
	  const ColumnVector D1 = get_global_value ("D1") . column_vector_value ();
	  const ColumnVector D3 = get_global_value ("D3") . column_vector_value ();
	  const ColumnVector B1 = get_global_value ("B1") . column_vector_value ();
	  const ColumnVector B3 = get_global_value ("B3") . column_vector_value ();
	  const ColumnVector r10 = get_global_value ("r10") . column_vector_value ();
	  const ColumnVector r30 = get_global_value ("r30") . column_vector_value ();
	  const ColumnVector R_1 = get_global_value ("R_1") . column_vector_value ();
	  const ColumnVector C = get_global_value ("C") . column_vector_value ();
	  const ColumnVector Si = get_global_value ("Si") . column_vector_value ();
	  const ColumnVector A = get_global_value ("A") . column_vector_value ();

	  if (! error_state)
	    {
	      double r3 = r1 + r2;

	      ColumnVector r (3);
	      r(0) = r1;
	      r(1) = r2;
	      r(2) = r3;

	      ColumnVector E1 (3, 0.0);
	      Matrix E3 (3, 2, 0.0);
	      ColumnVector Q (3, 0.0);
	      ColumnVector J (3, 0.0);

	      for (int i = 0; i < 3; i++)
		{
		  E1(i) = D1(i)*(exp(-2*B1(i)*(r(i)-r10(i)))-2*exp(-B1(i)*(r(i)-r10(i))));
		  E3(i,0) = D3(i)*(exp(-2*B3(i)*(r(i)-r30(i)))+2*exp(-B3(i)*(r(i)-r30(i))));
		  E3(i,1) = C(i)*(r(i)+A(i))*exp(-Si(i)*r(i));

		  int j = (r(i) > R_1(i));

		  Q(i) = (E1(i) + E3(i,j)) / 2;
		  J(i) = (E1(i) - E3(i,j)) / 2;
		}

	      retval = Q(0)+Q(1)+Q(2)-sqrt(0.5*((J(0)-J(1))*(J(0)-J(1))+(J(1)-J(2))*(J(1)-J(2))+(J(0)-J(2))*(J(0)-J(2))));
	    }
	  else
	    error ("Vau: error getting global values");
	}
      else
	error ("Vau: invalid arguments");
    }
  else
    print_usage ("Vau");

  return retval;
}