<font size="2"><span style="font-family: Arial;">I got my program to run. But

the results doesn't make any sence.<br><br>At t

=&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; 0&nbsp; y

=1.000000e+00&nbsp; 0.000000e+00&nbsp; 0.000000e+00, <br>At t

=&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; 0.01&nbsp; y =9.996000e-01&nbsp;

4.000000e-04&nbsp; 0.000000e+00, <br>At t =&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;

0.02&nbsp; y =9.992002e-01&nbsp; -4.720016e-02&nbsp; 4.800000e-02, <br>At t

=&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; 0.03&nbsp; y =7.722397e-01&nbsp;

-6.681768e+02&nbsp; 6.684045e+02, <br>At t

=&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; 0.04&nbsp; y =-4.466124e+07&nbsp;

-1.338934e+11&nbsp; 1.339381e+11, <br>At t

=&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; 0.05&nbsp; y =-1.793342e+24&nbsp;

-5.376439e+27&nbsp; 5.378233e+27, <br>At t

=&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; 0.06&nbsp; y =-2.891574e+57&nbsp;

-8.668938e+60&nbsp; 8.671830e+60, <br>At t

=&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; 0.07&nbsp; y =-7.517556e+123&nbsp;

-2.253763e+127&nbsp; 2.254515e+127, <br>At t

=&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; 0.08&nbsp; y =-5.081142e+256&nbsp;

-1.523326e+260&nbsp; 1.523834e+260, <br>At t

=&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; 0.09&nbsp; y =-inf&nbsp; nan&nbsp; inf,

<br>At t =&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; 0.1&nbsp; y =nan&nbsp;

nan&nbsp; nan, <br>At t =&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; 0.11&nbsp; y

=nan&nbsp; nan&nbsp; nan, <br>At t =&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;

0.12&nbsp; y =nan&nbsp; nan&nbsp; nan, <br>At t

=&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; 0.13&nbsp; y =nan&nbsp; nan&nbsp; nan,

<br>At t =&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; 0.14&nbsp; y =nan&nbsp;

nan&nbsp; nan, <br><br>Would you please help me check what's wrong with my

code?<br><br>#include "petscts.h"<br><br>/*The problem is from<br>&nbsp;*

chemical kinetics, and consists of the following three rate<br>&nbsp;*

equations:&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;

<br>&nbsp;*&nbsp;&nbsp;&nbsp; dy1/dt = -.04*y1 +

1.e4*y2*y3<br>&nbsp;*&nbsp;&nbsp;&nbsp; dy2/dt = .04*y1 - 1.e4*y2*y3 -

3.e7*(y2)^2<br>&nbsp;*&nbsp;&nbsp;&nbsp; dy3/dt = 3.e7*(y2)^2<br>&nbsp;* on the

interval from t = 0.0 to t = 4.e10, with initial<br>&nbsp;* conditions: y1 =

1.0, y2 = y3 = 0.*/<br><br><br>typedef struct {<br>&nbsp;&nbsp;&nbsp; PetscInt

m;<br>}AppCtx;<br><br>extern PetscErrorCode&nbsp;

FormJacobian(TS,PetscReal,Vec,Mat*,Mat*,MatStructure*,void*),Monitor(TS

ts,PetscInt step,PetscReal time,Vec u,void*ctx),<br>&nbsp;&nbsp;&nbsp;&nbsp;

FormFunction(TS,PetscReal,Vec,Vec,void*);<br><br>int main(int argc,char

**argv)<br>{<br>&nbsp;&nbsp;&nbsp; TS ts;<br>&nbsp;&nbsp;&nbsp; Vec

u;<br>&nbsp;&nbsp;&nbsp; PetscScalar&nbsp;&nbsp;&nbsp;

*u_localptr;<br>&nbsp;&nbsp;&nbsp; Mat J;<br>&nbsp;&nbsp;&nbsp; AppCtx

user;<br>&nbsp;&nbsp;&nbsp; PetscInt its,m;<br>&nbsp;&nbsp;&nbsp; PetscReal

dt,ftime;<br>&nbsp;&nbsp;&nbsp; PetscErrorCode ierr;<br>&nbsp;&nbsp;&nbsp;

PetscInitialize(&amp;argc,&amp;argv,PETSC_NULL,PETSC_NULL); &nbsp;&nbsp;&nbsp; m=3; &nbsp;&nbsp;&nbsp; dt=0.01; &nbsp; ierr = PetscOptionsGetInt(PETSC_NULL,"-m",&amp;m,PETSC_NULL);CHKERRQ(ierr); &nbsp; user.m=m; &nbsp; ierr=VecCreateSeq(PETSC_COMM_SELF,m,&amp;u);CHKERRQ(ierr); &nbsp; ierr=VecGetArray(u,&amp;u_localptr);CHKERRQ(ierr); &nbsp; u_localptr[0]=1.0; &nbsp; u_localptr[1]=0.0; &nbsp; u_localptr[2]=0.0; &nbsp; ierr = VecRestoreArray(u,&amp;u_localptr);CHKERRQ(ierr);&nbsp;&nbsp;&nbsp;  &nbsp; ierr = TSCreate(PETSC_COMM_WORLD,&amp;ts);CHKERRQ(ierr); &nbsp; ierr = TSSetProblemType(ts,TS_NONLINEAR);CHKERRQ(ierr); &nbsp; ierr = TSMonitorSet(ts,Monitor,&amp;user,PETSC_NULL); &nbsp; ierr = TSSetSolution(ts,u);CHKERRQ(ierr); &nbsp; ierr = TSSetRHSFunction(ts,FormFunction,&amp;user);CHKERRQ(ierr); &nbsp; ierr = MatCreate(PETSC_COMM_SELF,&amp;J);CHKERRQ(ierr); &nbsp; ierr = MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,m,m);CHKERRQ(i!

err); &nbsp; ierr = MatSetFromOptions(J); &nbsp; ierr = TSSetRHSJacobian(ts,J,J,FormJacobian,&amp;user);CHKERRQ(ierr); &nbsp; ierr = TSSetInitialTimeStep(ts,0.0,dt);CHKERRQ(ierr); &nbsp; ierr = TSSetDuration(ts,1000,8.e-1); &nbsp; ierr = TSSetFromOptions(ts);CHKERRQ(ierr); &nbsp; ierr = TSSetUp(ts);CHKERRQ(ierr); &nbsp; ierr = TSStep(ts,&amp;its,&amp;ftime);CHKERRQ(ierr); &nbsp;  &nbsp; printf("Number of timesteps = %d final time %4.2e\n",(int)its,ftime); &nbsp; ierr = VecDestroy(u);CHKERRQ(ierr); &nbsp; ierr = MatDestroy(J);CHKERRQ(ierr); &nbsp; ierr = TSDestroy(ts);CHKERRQ(ierr); &nbsp; ierr = PetscFinalize();CHKERRQ(ierr); &nbsp; return 0; } PetscErrorCode FormFunction(TS ts,PetscReal t,Vec X,Vec F,void *ptr) {  &nbsp; PetscErrorCode ierr; &nbsp; PetscScalar *x,*f; &nbsp; PetscInt n; &nbsp; n=3; &nbsp; ierr = VecGetArray(X,&amp;x);CHKERRQ(ierr); &nbsp; ierr = VecGetArray(F,&amp;f);CH!

 KERRQ(ierr);<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; f[0!

 ]=-0.04*

x[0]+1.0e4*x[1]*x[2]; &nbsp;&nbsp;&nbsp; f[1]=0.04*x[0]-1.0e4*x[1]*x[2]-3*1.0e7*x[1]*x[1]; &nbsp;&nbsp;&nbsp; f[2]=3*1.0e7*x[1]*x[1]; &nbsp;&nbsp;&nbsp; ierr = VecRestoreArray(X,&amp;x);CHKERRQ(ierr); &nbsp; ierr = VecRestoreArray(F,&amp;f);CHKERRQ(ierr); &nbsp; return 0; &nbsp;  } PetscErrorCode FormJacobian(TS ts,PetscReal t,Vec X,Mat *J,Mat *B,MatStructure *flag,void *ptr) { &nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; Mat jac=*J; &nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; PetscScalar v[3],*x; &nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; PetscInt row,col; &nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; PetscErrorCode ierr; &nbsp;&nbsp;&nbsp; &nbsp;&nbsp; ierr = VecGetArray(X,&amp;x);CHKERRQ(ierr); &nbsp;&nbsp;&nbsp; &nbsp;&nbsp; v[0]=-0.04; &nbsp;&nbsp;&nbsp; &nbsp;&nbsp; v[1]=0.04; &nbsp;&nbsp;&nbsp; &nbsp;&nbsp; v[2]=0.0; &nbsp;&nbsp;&nbsp; &nbsp;&nbsp; row=0; &nbsp;&nbsp;&nbsp; &nbsp;&nbsp; col=0; &nbsp;&nbsp;&nbsp; &nbsp;&nbsp; ierr = Ma!

tSetValues(jac,3,&amp;row,1,&amp;col,v,INSERT_VALUES);CHKERRQ(ierr); &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; v[0]=1.0e4*x[2]; &nbsp;&nbsp;&nbsp; &nbsp;&nbsp; v[1]=-1.0e4*x[2]-6.0e7*x[1]; &nbsp;&nbsp;&nbsp; &nbsp;&nbsp; v[2]=6.0e7*x[1]; &nbsp;&nbsp;&nbsp; &nbsp;&nbsp; col=1; &nbsp;&nbsp;&nbsp; &nbsp;&nbsp; ierr = MatSetValues(jac,3,&amp;row,1,&amp;col,v,INSERT_VALUES);CHKERRQ(ierr); &nbsp;&nbsp;&nbsp; &nbsp;&nbsp; v[0]=1.0e4*x[1]; &nbsp;&nbsp;&nbsp; &nbsp;&nbsp; v[1]=-1.0e4*x[1]; &nbsp;&nbsp;&nbsp; &nbsp;&nbsp; v[2]=0.0; &nbsp;&nbsp;&nbsp; &nbsp;&nbsp; col=2; &nbsp;&nbsp;&nbsp; &nbsp;&nbsp; ierr = MatSetValues(jac,3,&amp;row,1,&amp;col,v,INSERT_VALUES);CHKERRQ(ierr); &nbsp; ierr = MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); &nbsp; ierr = VecRestoreArray(X,&amp;x);CHKERRQ(ierr); &nbsp; ierr = MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); &nbsp; return 0; } PetscErrorCode Monitor(TS ts!

,PetscInt step,PetscReal time,Vec u, void *ctx) { &nbsp;!

  PetscSc

alar *y;<br>&nbsp; <br>&nbsp; VecGetArray(u,&amp;y);<br>&nbsp; PetscPrintf(PETSC_COMM_WORLD,"At t =%11g&nbsp; y =%e&nbsp; %e&nbsp; %e, \n",time,y[0],y[1],y[2]);<br>&nbsp; VecRestoreArray(u,&amp;y);<br>&nbsp; return 0;<br>}<br><br><br></span></font>