#
# Cross diffusion epidemiology model
#
# Delete previous workspaces
  rm(list=ls(all=TRUE))
#
# Access ODE integrator
  library("deSolve");
#
# Access functions for numerical solutions
  setwd("f:/comp3/wiley_download/chap7");
  source("pde_1.R");
  source("dss004.R");
  source("dss044.R");
#
# Level of output
#
#   ip = 1 - graphical (plotted) solutions 
#            (u1(x,t), u2(x,t)) only
#
#   ip = 2 - numerical and graphical solutions
#
  ip=2;
#
# Parameters
  a=0.05;     b=0.006; 
  alpha=0.06; beta=0.0056; gamma=0.04; 
  m=2;        n=1; 
  d11=100;    d12=6;       d22=3;
#
# Revised values
  d11=3; a=0.006; alpha=0.006; gamma=0.006; 
#
# Select case
  ncase=3;
#
# Initial conditions
  nx=41;
  u0=rep(0,2*nx);
#
# ncase=1
  if(ncase==1){
    for(i in 1:nx){
      u0[i]=1;
      u0[i+nx]=0;
    }
    d12=0;
  }
#
# ncase=2,3
  if(ncase>1){
    for(i in 1:nx){
      u0[i]=1;
      if(i<=5){
        u0[i+nx]=1;
      }else{
      u0[i+nx]=0; 
      }
    }
    if(ncase==2){d12=0;}
    if(ncase==3){d12=6;}
  }
  ncall=0;  
#
# Write selected parameters
  cat(sprintf("\n\n ncase = %2d  d11 = %4.2f  d12 = %4.2f  d22 = %4.2f",
              ncase,d11,d12,d22));
#
# Write heading
  if(ip==1){
    cat(sprintf("\n Graphical output only\n"));
  } 
#
# Grid (in x)
  xl=0;xu=20;
  x=seq(from=xl,to=xu,by=(xu-xl)/(nx-1));     
#
# Independent variable for ODE integration
  nout=5;t0=0;tf=20;
  tout=seq(from=t0,to=tf,by=(tf-t0)/(nout-1));
#
# ODE integration
  out=lsodes(y=u0,times=tout,func=pde_1,parms=NULL)
  nrow(out)
  ncol(out)
#
# Arrays for plotting numerical solution
  u1_plot=matrix(0,nrow=nx,ncol=nout);
  u2_plot=matrix(0,nrow=nx,ncol=nout);
  for(it in 1:nout){
    for(ix in 1:nx){
       u1_plot[ix,it]=out[it,ix+1];
       u2_plot[ix,it]=out[it,ix+1+nx];
    }
  }
#
# Display numerical solution
  if(ip==2){
    for(it in 1:nout){
      cat(sprintf(
      "\n    t       x     u1(x,t)     u2(x,t)\n"));
      for(ix in 1:nx){
        cat(sprintf("%5.1f%8.1f%12.3f%12.3f\n",
        tout[it],x[ix],u1_plot[ix,it],u2_plot[ix,it]));
      }
    }
  }
#
# Calls to ODE routine
  cat(sprintf("\n\n ncall = %5d\n\n",ncall));
#
# Plot u1
  par(mfrow=c(1,1));    
  matplot(x=x,y=u1_plot,type="l",xlab="x",
          ylab="u1(x,t), t=0,5,...,20",xlim=c(xl,xu),lty=1,
          main="u1(x,t); t=0,5,...,20;",lwd=2);
#
# Plot u2
  par(mfrow=c(1,1));    
  matplot(x=x,y=u2_plot,type="l",xlab="x",
          ylab="u2(x,t), t=0,5,...,20",xlim=c(xl,xu),lty=1,
          main="u2(x,t); t=0,5,...,20;",lwd=2);