Lecture 2.

 

First order differential equations in dilution problems and other situations

 

by Reinaldo Baretti Machνn

http://www1.uprh.edu/rbaretti

http://www1.uprh.edu/rbaretti/methodsoftheoreticalphysics.htm

http://www1.uprh.edu/rbaretti/methodsoftheoreticalphysicsPart2.htm

http://www1.uprh.edu/rbaretti/methodsoftheoreticalphysicsPart3.htm

 

 

Para preguntas y sugerencias escriba a :

reibaretti2004@yahoo.com

 

 

Free counter and web stats  

 

 

Reference :

1. Mathematical Techniques for Biology and Medicine (Dover Publications- Paperback)

by William Simon ; http://www.amazon.com/Mathematical-Techniques-Biology-Medicine-William/dp/0486652475/ref=sr_1_1?ie=UTF8&s=books&qid=1223731765&sr=8-1

 

2. http://geocities.com/serienumerica4/Lecture1Biolologymedicine.doc

 

 

Examples and problem from:

Chapter 2 . Reference 1. The problem numbers refer to the textbook listing.

 

Section 1. One compartment dilution problem.

 

Fig 1.

A bucket contains say V=10 liters of water. The volume is kept constant becaues the inflow F =2 liter/minute  is adjusted is equal to the outflow F at the bottom. Thres moles of a tracer Q(0)= 3 moles are injetct at t=0 and we assume it is instantly mixed throughout the volume V.

 

At any instant ,the concentration  C(t) (moles/liter) = Q(t) /V  .

 

The initial concentration C(0) = 3 moles/10 liters =0.3 moles/liter

 

Establish a first order differential equation and solved it.

 

In an interval  ∆t (minutes) the number of moles diminishes by ∆Q .

 

A finite difference equation is

 

 

∆Q = - ∆t  F C(t)                        (1)

moles ~ minutes * (liters/minute) * (moles/liter)

We can write an equation for either Q or C.

 

Substitute   ∆Q = V ∆C  above an get

∆C  = - ∆t ( F/V) C(t)  . The differential equation follows from invoking the

 

limit (∆t→0) (∆C / ∆t ) = - (F/V) C

 

  i.e       dC/dt = -(F/V) C                                          (2)

 

This is of the form dy/dt = -k y (see chapter 1) and the solution is

 

       y(t) = y(0)exp(-kt)  . Hence the solution is,

 

C(t) = C(0) exp( -(F/V) t)   . The numerical values are

C(0)= 0.3 moles/liter   , F/V =k or λ = 2(liters/min)/10(liters) = .2 /minute

Tscale = 1/λ = 5 minutes The plot of C(t) is shown next.

 

 

 

Plots of C(t) analytical and the numerical result from an integration.

 

 

The differential equation

dC/dt = -(F/V) C     leads to the finite difference equation

 

(Cn – Cn-1 )/∆t = -(F/V) Cn-1  wrere Cn = C( n ∆t)

 

Solving for Cn

 

Cn =Cn-1  - ∆t (F/V) Cn-1                                            (3)

 C0 = C(t=0)   is given  , then the solution is constructed from the previous

values of C which are known.

 

C1 =C0  - ∆t (F/V) C0 ,

C2 =C1  - ∆t (F/V) C1 …

 

Now  ∆t has to be small relative to the time scale 1/λ .

The following MATLAB code integrates numerically the DE.

The value of ∆t is chosen to be small in comparison with Tscale.

C analytical  i.e C(0) exp(-λt) and C numerical are plotted together. They are practically identical within the resolution of the graph.

 

tscale = 5 minutes    , dt= 2.0 E-2 minutes

 

MATLAB code

%numerical solution of dC/dt= -(F/V)*C(t)

 

F=2 ; V=10 ; Q0= 3; c(1)=Q0/V   ;

lambda=F/V; tscale=1/lambda  ,

cexact=c(1)*exp(-lambda*t);

ti=0 ; tf = 4*tscale ; nstep=1000;

dt=(tf-ti)/nstep ,

% k is the index kounter it starts with k=1 at t=0

% c(1) is given as initial condition , we wil find c*2) ,

% c(3)...etc

k=0 ;

for x=[ti+dt:dt:tf]

k=k+1; 

c(k+1)=c(k) +dt*(-(F/V)*c(k));

 

end

 

t=[ti:dt:tf];

plot(t,c,t,cexact), xlabel('t(minutes)'),ylabel('C,cexact')

 

 

Other examples of exponential decay.

 

Section 2.

 

Consider a metabolite of Q0 moles to which a small quantity of radioactively

 

labeled q(0) is added at t=0. ( i.e q(0) << Q0 ) . The total quantity is

 

               

                         Qtotal = Q0 + q(t) ≈Q0 .

q decays at a fixed rate  of   R ~ moles/time

 

Before writing an equation for ∆q  let’s rewrite eq (1)

 

∆Q = - ∆t  F C(t)                        as

 

change in q =(moles)= - ∆t (time) R (moles/time) (radioactive moles/total moles)

 

∆q = - ∆t R q(t)/Q0  from which we can pass to the DE

 

dq/dt = -k q      k= R/Q0 ~ 1/time                           (4)

 

Radioactive decay

It is empirically found that from a total of N atoms, there is a decay ,  ∆N<0 , in an interval ∆t such that

 

         ∆N =- λ N ∆t , λ~ 1/time      .  The corresponding DE is

 

  dN/dt = - λ N(t)    . Let the initial condition be N(t=0) = N0 , then

    N = N0 exp( - λt)

 

Absorption of  X- rays

The intensity I ~ energy/(time-area)  , of X-rays after penetrating a material slab of thickness ∆x  is diminished acording to

∆I = - D I(x) ∆x    , D ~ 1/length

 

this leads to the DE   dI/dx = - D I  . Let the initital condition

be I(x=0) = I0 , then    I (x) = Io exp(-Dx)  

 

Sec 5. One compartment with constant injection   (see figure)

Assume the intial concentration of substance Q is zero , C(0)=0.

 

A volume V~ liter is kept fixed by a constant inflow flow and outflow of F ~

liters /time . The quantity of Q(t) moles of a  tracer suffers in an interval

a loss 

 

 ∆Qloss = - ∆t F Q/V ~ - ∆t ( volume/time) * (moles/volume).

 

At the same time a quantity  ∆t R ~ time *(moles/time) is injected instantly.

 

R is the rate of injection in moles/time .

 

∆Qgain = ∆t R   . Adding the loss and gain parts of Q ,

 

 

∆Q = ∆t (  R  -  F Q/V) =  ∆t (  R  -  F C(t) ).

 

Dividing by  ∆t V and taking the limit ∆t→0 , results in the DE

 

d C /dt = ( R/V – (F/V) C(t) )                 ,              (5)

 

with initial condition C(t=0) =0.

 

Remember R , V , F are parameters. They assume  constant values in a given example. The independent variable is time , the dependent variable is the concentration C(t) ~ moles/volume.

 

The DE equation (5) is of a type called separable. One can separate the variables  C and t , in two different integrals.

 

∫ dC /(b -a C(t) )       =  ∫ dt  , where , b=R/V   a= - F/V .

 

If you go to a table of integrals you will find the indefinite integral in a form  like this ;

 

∫ dx /(b+ax) = (1/a) ln (b+ax)   . Hence

 

∫ dC / { (R/V) – (F/V) C(t) }    =

      -(1/(F/V) ) ln [ (R/V) – (F/V) C(t)]  .

 

Equating this to the right hand side integral  gives

 

the indefinite integral

 

-(1/(F/V) ) ln [ (R/V) – (F/V) C(t)]  = t + A

 

A is the integration constant that can be related to the initial condition as we shall show.

 

It is convenient to rewrite the integral as

 

 ln [ (R/V) – (F/V) C(t)]  = - (F/V) t  + A*    where A* = (F/V) A is still a constant. To simplify matters call A all constants that originate with the original A.

 

Then

      (R/V) – (F/V) C(t) = A exp(-(F/V)t )   ,  (another A)

and solving for C(t)

 

C(t) = (V/F) { (R/V) - A exp(-(F/V)t )  }                       (6).

 

Apply now the initial condition to determine A

 

C(0)= 0 = (V/F) { (R/V) – A }   , thus A =R/V  and finally

 

C(t) = (R/F) { 1- exp( -(F/V) t)  }                                   (7)

which dimensionally checks ,

 

moles/volume = ( moles/time)( time/volume) .

At t=” infinty”     C(∞) = R/F = constant  .

The time scale of the problem is T scale  = V/F

 

Example : Plot of eq (7) . Such plot is called an inverted exponential.

 set   % R=.1 mol sec , F =2 liter/(60 sec) , V=10 liters

 

 

 

Matlab code

% One compartment with constant injection

% R=.1 mol sec , F =2 liter/(60 sec) , V=10 liters

R=0.1 , F =2/60 , V=10 ,

tscale=V/F; tfinal=3*tscale ;dt = tfinal/100;

t=[0:dt:tfinal];

c=(R/F)*(1-exp(-t/tscale));

plot(t,c) , xlabel('t(sec)') ,...

    ylabel('C(t)~moles/liter ' )

 

The two compartment series dilution

The one compartment dilution problem have already been solved

and the solution is C1(t)=C1(0) exp( -(F/V1) t)   where the subscript 1 refers to compartment #1. The initial condition was C1(0)= Q0 /V1  i.e a total amount of tracer Q0 is injected and mixed instantly at t=0.

Volumes V1 and V2 are kept constant by insuring that F = F1 = F2 .

 

The problem is to set up a differential equation for C2(t) and to find its solution. Remember C2(t) = Q2(t) /V2 and Q2 ~ moles in the second compartment. The   initial condition is C2(t) =0.

 

The quantity  ∆Q2 ~ moles ,in compartment 2 within an interval ∆t is increased by the amount 

 

∆Q2 gain = inflow rate * C1 * ∆t~ (liter/sec)(moles/liter) (sec)

             = F  C1 ∆t

 

The loss is  ( negative quantity) 

∆Q2 loss = - outflow rate * C2 * ∆t~ (liter/sec)(moles/liter) (sec).

 

             =  - F C2 (t) ∆t

 

So the net change is

 

∆Q2 = F  C1 ∆t - F C2 (t) ∆t  and the corresponding DE is

 

d Q2 /dt = F  C1  - F C2 (t)  divide by V2 to obtain the DE for C2 ,

 

dC2/dt = (F/V2 ) ( C1(t) – C2 (t) )               .  (8) 

 

Rearranging (8) , and introducing the solution for C1  get

 

dC2/dt  +( F/V2 )  C2 (t) = (F/V2 ) C1(0) exp( -(F/V1) t )   . (9)

 

This is a first order linear DE  of the form

 

   dy/dx +  P(x) y  = Q(x)                                               (10)

 

which has the general solution

 

y = exp(-I)  ∫ Q(x) exp(I) dx +  B exp(-I)                         (11)

 

B is an integration constant  and  I = ∫ P(x) dx  .   

 

To avoid the burden of too many symbols let’s put some data in (9).

 

Set , C1 (0)=0.3moles/liter   , V1= 10 liters  , V2= 5 liters ,

 

F= 2liters/min= 2liters/(60sec) = 3.33E-2 liters/sec.

There are two time scales ;  (F/V1) = 3.33E-3 sec-1 ; Tscale1 = V1/F= 300 sec   and (F/V2) = 6.66E-3 sec-1 ; Tscale2 = V2/F = 150 sec.

 

 With the proposed data and the value of  C1(0)=.3 moles/liter  , eq (9) becomes

 

dC2/dt  + 6.67E-3 C2 (t) = 6.67E-3 *(.3) * exp( -3.33E-3 t )  . (12)

 

Identify P(t) = 6.67E-3 ~ 1/sec  and

Q(t) = 2.00E-3 *exp(-3.33E-3 t)   ~    moles/(liter-sec)

 

Proceed to implement the method of eq  (11)

I = ∫ P(t) dt  = 6.67E-3 * t                             (13)

 

∫Q(x) exp(I) dx = ∫ 2.00E-3 *exp(-3.33E-3 t)*exp(6.67E-3  t) dt   

       

           =  ∫  2.00E-3 *exp(+3.33E-3 t) dt

            = ( 2.00E-3/3.33E-3)  exp( +3.33E-3 t) 

 

 

            =  6.01E-1 exp( 3.33E-3 t) .                                            (14)

 

Inserting (13), (14) in (11) results in

 

C2(t) = exp(-6.67E-3 t) * 6.01E-1 exp( 3.33E-3 t) + Bexp(-6.67E-3 t)

 

        = 6.01E-1 exp(-3.33E-3t) + B exp(-6.67E-3 t)  .

Now invoke the initail condition

 

C2(0)= 0 = 6.01E-1 + B   ; thus B= -6.01E-1  and

 

C2(t) = 6.01E-1 { exp(-3.33E-3t) – exp (-6.67E-3 t)  }                 .  (15)

 

Applying (11) but keeping the data (the parameters) unspecified gives,

 

C2(t) = [ V1 C1(0)/(V1-V2 )] {exp(-(F/V1)t ) – exp(-(F/V2)t )   }  . (16)

 

         

 

 

 

Plot of C1 , C2

 

Matlab code

 

% two compartment dilution

tscale=1/6.67E-3 ; tf=5*tscale;

dt=tf/200;

t=[0:dt:tf];

c1=.3*exp(-3.33E-3*t);

c2=6.01E-1*(exp(-3.33E-3*t)-exp(-6.67E-3*t) );

plot(t,c2,t,c1), xlabel('t~sec') ,...

    ylabel('C1 ,C2~moles/liter')

 

 

Roughly speaking without solving the DE we can anticipate  the following

 

results. The  3 moles are emptied from compartment one in roughly 300 seconds ( tscale 1). If there is no loss in compartment #2 the concentration would be 3 moles/5 liters = 0.6 moles/liter  but in 300 sec there are two half lifes of tscale2 (tscale2=150 sec) . Thus the concentratio reduces to (1/2)2=1/4 . So at t=300 sec the concentration at compartment # 2 ,is at a maximum and roughly equals   (1/4)(.6moles/liter) = 0.15 moles/liter .

 

A  check of the plot above bears out this estimate.

 

Numerical integration of eq.(8).

 

We expand next on the procedure already used ( in these notes) to integrate numerically a DE by finite differences.

 

The equation

 

dC2/dt  +( F/V2 )  C2 (t) = (F/V2 ) C1(0) exp( -(F/V1) t )        (17)

has initial condition C2(t=0) = 0 . Suppose you want the solution in intervals

 

of size ∆t where ,  ∆t = ( tfinal – tinitial) /nsteps , and ∆t << tscale1 or

∆t << tscale2 , whichever is the samller one. 

 

 

 

 The finite difference generic  solution is ,

 

C2 (tn) = C2 (tn-1) +

                    ∆t  {- ( F/V2 )  C2 (tn-1) + (F/V2 ) C1(0) exp( -(F/V1) tn-1 ) } . (18)

 

note that every term  in the right hand is known.

 

The Fortran given below uses this method. The plot is given next , which , of course is identical to the previous one.

 

Same initial data C1(0)= 0.3 moles/liter , c2(0)=0,  but V1=V2=10 liters.

 

Fortran code (two compartment series dilution)

c the two compartment series dilution

c F~liters/min ,V ~liters

      data c10 /.3 /

      data F, V1,V2/2. ,10.,  5./

      c1(t) = c10*exp(-F*t/v1)

c  f is changed to liters /second by the next declaration

      f=f*(1./60.)

      tscale1=V1/f

      tscale2=V2/f

      tsmall=amin1(tscale1,tscale2)

      tlarge=amax1(tscale1,tscale2)

      dt= tsmall/50.

      nstep=int(3.*tlarge/dt)

      c20=0.

      kp=int(float(nstep)/70.)

      kount=kp

      print 100 , 0. , c1(0.) ,c20

      do 10 i=1,nstep

      t=dt*float(i)

      c21= c20 + dt*(-(F/v2)*c20 +(f/v2)*c1(t-dt) )

      if(i.eq.kount)then

      print 100 ,t ,c1(t),c21

      kount=kount+kp

      endif

      c20=c21

10    continue

100   format(1x,'t(sec),c1, c2 ~moles/liter=',3(3x,e10.3))

      stop

      end

     

 

Section 8. Mother Daughter Radioactive Decay

 

This problem is analogous to the two compartment dilution series.

 

A mother (M) nuclei species decays like

 

dM/dt =  - λ1 M                                                        ( 19  ).

 

The daughter nuclei ( D) are born at the rate  + λ1 M    and decay at a rate

- λ2 D  , therefore the DE for the daughter nucei is

 

d D/dt =   λ1 M(t)     - λ2 D (t)                                  (20) 

The equation is analogous ( though not the same) to

 

dC2/dt = (F/V2 ) ( C1(t) – C2 (t) )      if we make C1=M , C2= D

 

and    λ1 = λ2 = (F/V2) .    But in general  λ1 ≠ λ2  .

 

We need two initial conditions to solve (19) and (20) , set M(0)= 1mol

D(0)=0  and  tscale1 =1/ λ1 = 5 minutes , tscale2 = 1/λ2 =10 minutes.

 

We can modify slightly the above FORTRAN code to integrate , both equations simultaneously. The finite difference solutions to (19)-(20) are

 

M(tn ) =  M(tn-1) - ∆t λ1 M(t n-1)                            (21)

 

                                           

D( tn )  =  D(tn-1)  + ∆t (  λ1 M(tn-1 )     - λ2 D (tn-1)  )    . (22)

 

Solution of (21)-(22) with the above parameters

 

FORTRAN code

c Fortran code for mother daughter radioactive decay

      real m0 ,m1 ,lambda1, lambda2

      data m0 ,D0 /1.,0. /

c time scales are in minutes

      data tscale1 , tscale2/5.,10. /

      lambda1=1./tscale1

      lambda2=1./tscale2

      tsmall=amin1(tscale1,tscale2)

      tlarge=amax1(tscale1,tscale2)

      dt= tsmall/50.

      nstep=int(3.*tlarge/dt)

      kp=int(float(nstep)/70.)

      kount=kp

      print 100 , 0. , m0 ,d0

      do 10 i=1,nstep

      t=dt*float(i)

      M1 =  M0 - lambda1*dt*M0

      D1=d0+ dt*(lambda1*m0-lambda2*d0)

      if(i.eq.kount)then

      print 100 ,t ,m1,d1

      kount=kount+kp

      endif

      m0=m1

      d0=d1

10    continue

100   format(1x,'t(min),mother,daughter~moles=',3(3x,e10.3))

      stop

      end

 

 

    

 

 

Section 9. Circular reactions

 

                 

The concentration A has a gain in an interval ∆t   of

 

∆t   k3 C C ‘   (third line)   and a loss of   - ∆t  k1 A A ‘   (first line). The balance is

 

∆A = ∆t   k3 C C ‘   - ∆t   k1 A A ‘   , which leads to the DE

 

dA/dt = k3 C C ‘   -  k1 A A’                                       (23)

 

 

the ki ~ ( 1/time)*(1/concentration) .

Similarly

 

dB/dt = k1 A A ‘   -  k2 B B’                                     (24)

 

dC/dt = k2 B B ‘   -  k3 C C’   .                                 (25)

 

To simplify matters suppose A’ , B’ and C’ are practically constants and can be absorbed into new constants , capital K1 , K2 , K3 .

 

The capital K’s are ;  K1 ≡ k1 A ‘    , K2 ≡ k2 B ‘ ,   K3 ≡ k3 C ‘ ~1/time  . The

equations take the new form

 

 

dA/dt = K3 C   -  K1 A   ,                                       26-a

 

dB/dt= K1A   -K2 B  ,                                             26-b

 

dC/dt = K2 B – K3 C   .                                         26-c

 

We will solve numerically this system of three first order coupled DE.

 

  Let  tscale1= 10 minutes =1/K1  , tscale2=5minutes =1/K2 ,

 

Tscale3= 3minutes=1/K3  .  A(0)=1.0 moles , B(0)=0 moles , C(0)=0 moles

 

The finite difference solutions to eqs 26 a to  26 c  are

 

A( tn ) =A(tn-1) +  ∆t { K3 C (tn-1)  -  K1 A(tn-1 )  }                           

B( tn ) =B(tn-1) +  ∆t { K1 A (tn-1)  -  K2 B(tn-1 )  }   

 C( tn ) =C(tn-1) +  ∆t { K2 B (tn-1)  -  K3 C(tn-1 )  }     

and ∆t has to be small relative to the smallest of the Tscale. In the present

case     ∆t << Tscale 3.   

 

The next FORTRAN code expands the mother –daughter code to deal with the circular reaction.  The results are plotted next.

 

We see that in about 10 minutes a steady state is reached with constant concentrations

A=.56   , B= .28   , C= .16 .        

                      

 

Problem 3. A substance Q~ moles , diffuses from the compartment of volume V at a rate proportional to the product  of the diffusion constant K ~ volume/time and the difference in concentration with the environment (C- Cenvironment) ~ moles/volume. The initial condition is C(0)= Q0/V.

Also assume C > Cenvironment hence there is a loss of Q in the compartment,

V is kept constant , the mechanism is not shown in the drawing.

Then in an interval  ∆t , a loss occurs

 

∆Q = - ∆t K (C- Cenvironment) ~ time( volume/time)*(moles/volume)

 

Dividing by the volume V and by ∆t we obtain the differential equation

 

          dC/dt =  -(K/V)  (C(t) - Cenvironment) .

 

Problem 4. Assume that in the compartement of Prob. 3 a quantity ∆Qin

is pumped given by flow rate R(moles/time) * ∆t .

 

The net ∆Q is now

 

∆Q  =  - ∆t K (C- Cenvironment) + R ∆t  and the DE is

 

dC/dt = - (K/V)(C- CE ) +  R/V

 

Example : Solve the DE dC/dt = - (K/V)(C- CE ) +  R/V .

   Let K= F=2liter/minute   , V=10 liter  ,C(0)=0.3 moles/liter ,

 

Cenvironmnet=.05 moles/liter   , R =.1mole/minute.

 

 

Tscale = 1/(F/V) =10/2 = 5 minutes .

 

The finite difference solution is 

 

C ( tn) = C(tn-1)  + ∆t { - (K/V)(C(tn-1) - CE )  + R/V }.

FORTRAN code

      real k

      data  c0 , cE/.3, .10/

      data k,v, R/2.,10.,0.1/

c time scales are in minutes

      tscale=V/K

      dt= tscale/100.

      nstep=int(6.5*tscale/dt)

      kp=int(float(nstep)/70.)

      kount=kp

      print 100 , 0. , c0

      do 10 i=1,nstep

      t=dt*float(i)

      c1 =  c0 +dt*( -(K/v)*(c0-ce) +R/V)

      if(i.eq.kount)then

      print 100 ,t ,c1

      kount=kount+kp

      endif

      c0=c1

10    continue

100   format(1x,'t(min),C~ moles/liter=',2(3x,e10.3))

      stop

      end

 

 

 

Problem 5 . In a one compartment dilution the concentration falls with  time according to the data :

 t(sec)

C (moles/liter)

0

.024

1

.011

2

4.8E-3

3

2.4E-3

4

1.0E-3

 

If Q(0)=0.1 mole , find F ~liters/sec   and  V ~liters

 

C(0) =.024 ≡ Q(0)/V  , hence V= .1/.024 = 4.17 liters  

 

If we assume a perfect fit , which may not be the case ,

 

then  C(t) = C(0) exp(- (F/4.17) t) . Any value from the table will serve to solve for F. At t=1 sec

 

.011 = .024 exp( -.240F)  . Taking the natural log

 

ln(.011/.024) = -.780 = -.240F

 

F = 3.25 liters/second

 

Problem 6.  

In a radioactive  process half the  atoms decayed in one day.What fraction existed at 10 hours. We know that the half life (t1/2 ), is related to the rate λ by

 

                      λ = .693 / t1/2 = .693 /24 hours = 2.89E-2 hour-1 .

 

After 10 hours the fraction of atoms left is

 

N/N0 = exp( - λ t)  = exp( -2.89E-2 (10) ) = .749  .

 

 

Problem 8.  Substance A is created at a rate R ~moles/sec and destroyed at the rate   - K(1/sec)  * A(moles)  ~ moles/sec. Let  the initial condition be A(0) =0.

In an interval ∆t  the net change in A is

 

∆A =  R ∆t – K A ∆t

The corresponding DE is

 

dA/dt = R – K A  .

 

A steady state is reached when dA/dt =0 . This gives A(t=∞) = R/K.

Again the solution is an inverted exponential as mentioned in a previous example of the one compartment dilution with constant injection.

 

Data:  R= .5 moles/sec , K = 0.1/sec , A(0)=0.

The asymptotic value of A is  A (t=∞) = R/K=5 moles

 

MATLAB code

 

% R~moles/sec K ~(moles/mole-sec)

A(1)=0; R =.5 ; K= .1   ;

tscale= 1/K ;tfinal=3*tscale ;nstep=2000;

dt=tfinal/nstep;

n=0;

 t=[dt:dt:tfinal];

n=n+1;

A(n+1)=A(n) + dt*(R-K*A(n));

 

t=[0:dt:tfinal];

plot(t,A),xlabel('t~sec'),ylabel('A~moles')

 

 

 

 

 

 

END OF CHAPTER II