`help/text/DIFFALL`:=TEXT(
`HELP FOR:numerical solution of the diffusion partial differential`,
`         eqaution for any one of FTCSexp (foward time, centred `,
`         space, explict), FTCSimp or CRANK-NICHOLSON techniques`,
` `,
`CALLING SEQUENCE:`,
`   DIFFALL();`,
` `,
`SYNOPSIS:`,
`-no parameters are required as input`,
`-all input is prompted`,
` `,
`EXAMPLE:`,
`To solve the diffusion equation numerically using the FTCSimp`,
`scheme, select this scheme when prompted by typing`,
`         FTCSimp`,
`The numerical solution for the initial square wave evolving`,
`according to the diffusion equation on the specified mesh is then`,
`calculated and shown as a 3 dimensional plot`
):

#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
#>>>>>

DIFFALL:=proc(a,k,NN,JJ,SMAT,edgeS,edgeL,schm)
       local scm,VEC,GMAT,EMAT,MMAT,i,ii,ij;

#compare implicit, explicit FTCS and Crank-Nicholson
#for solving the DIFFUSION equation
#now entering the iterative scheme as a parameter (readstat)
#PLOTS and LINALG are required


#NN:=readstat(`no of time steps__`);
#JJ:=readstat(`no of space steps__`);
VEC:=array(1..JJ);
EMAT:=array(1..JJ,1..JJ);GMAT:=array(1..JJ,1..JJ);MMAT:=array(1..JJ,1..JJ);

#initialise matrices
for i from 1 to NN do
  for ii from 1 to JJ do
    SMAT[i,ii]:=0.0;
  od;
od;

for i from 1 to JJ do
  for ii from 1 to JJ do
    MMAT[i,ii]:=0.0;EMAT[i,ii]:=0.0;GMAT:=0.0;
  od;
od;

#starting wave
#edgeS:=readstat(`start of square wave__`);
#edgeL:=readstat(`end of square wave__`);

for i from 1 to JJ do VEC[i]:=0; od;
#define initial wave (square)
for i from edgeS to edgeL do VEC[i]:=1;SMAT[1,i]:=1; od;


#what kind of iterative scheme??????????
#EMAT and MMAT are to be extracted from the differencing schemes
#schm:=readstat(`differencing scheme? FTCSexp, FTCSimp, CRANK`);
scm:=op(1,schm)-op(2,schm);

#k:=readstat(`diffusion coefficient/function (non-linear)`);
#a:=readstat(`dt/dx**2__`);

#make MMAT and EMAT
for i from 2 to JJ-1 do

#3 substitutions for each
#1.diff(K..)=0 to remove unwanted terms
#2.K(u)=u^2 for example to insert diffusion function
#3.vec[]=u to insert current value of wavefunction

GMAT[i,i]:=-diff(scm,u[j,n]);
#1st sub
GMAT[i,i]:=subs(diff(K(u[j,n]),u[j,n])=0,GMAT[i,i]);
#2nd sub
GMAT[i,i]:=subs({K(u[j,n])=k(u[j,n]),K(u[j-1,n])=k(u[j-1,n]),
                 K(u[j+1,n])=k(u[j+1,n])},GMAT[i,i]);


GMAT[i,i+1]:=-diff(scm,u[j+1,n]);#print(1,");
GMAT[i,i+1]:=subs(diff(K(u[j+1,n]),u[j+1,n])=0,GMAT[i,i+1]);#print(2,");
GMAT[i,i+1]:=subs({K(u[j,n])=k(u[j,n]),K(u[j-1,n])=k(u[j-1,n]),
                 K(u[j+1,n])=k(u[j+1,n])},GMAT[i,i+1]);#print(3,");

GMAT[i,i-1]:=-diff(scm,u[j-1,n]);
GMAT[i,i-1]:=subs(diff(K(u[j-1,n]),u[j-1,n])=0,GMAT[i,i-1]);
GMAT[i,i-1]:=subs({K(u[j,n])=k(u[j,n]),K(u[j-1,n])=k(u[j-1,n]),
                 K(u[j+1,n])=k(u[j+1,n])},GMAT[i,i-1]);
#EMAT 2nd
EMAT[i,i]:=diff(scm,u[j,n+1]);
EMAT[i,i+1]:=diff(scm,u[j+1,n+1]);
EMAT[i,i-1]:=diff(scm,u[j-1,n+1]);
od;
#some extras!!!!!!!!!!!!!!!!!!!!
GMAT[1,1]:=-diff(scm,u[j,n]);
GMAT[1,1]:=subs(diff(K(u[j,n]),u[j,n])=0,GMAT[1,1]);
GMAT[1,1]:=subs({K(u[j,n])=k(u[j,n]),K(u[j-1,n])=k(u[j-1,n]),
                 K(u[j+1,n])=k(u[j+1,n])},GMAT[1,1]);

GMAT[1,2]:=-diff(scm,u[j+1,n]);
GMAT[1,2]:=subs(diff(K(u[j+1,n]),u[j+1,n])=0,GMAT[1,2]);
GMAT[1,2]:=subs({K(u[j,n])=k(u[j,n]),K(u[j-1,n])=k(u[j-1,n]),
                 K(u[j+1,n])=k(u[j+1,n])},GMAT[1,2]);

GMAT[JJ,JJ-1]:=-diff(scm,u[j-1,n]);
GMAT[JJ,JJ-1]:=subs(diff(K(u[j-1,n]),u[j-1,n])=0,GMAT[JJ,JJ-1]);
GMAT[JJ,JJ-1]:=subs({K(u[j,n])=k(u[j,n]),K(u[j-1,n])=k(u[j-1,n]),
                 K(u[j+1,n])=k(u[j+1,n])},GMAT[JJ,JJ-1]);

GMAT[JJ,JJ]:=-diff(scm,u[j,n]);
GMAT[JJ,JJ]:=subs(diff(K(u[j,n]),u[j,n])=0,GMAT[JJ,JJ]);
GMAT[JJ,JJ]:=subs({K(u[j,n])=k(u[j,n]),K(u[j-1,n])=k(u[j-1,n]),
                 K(u[j+1,n])=k(u[j+1,n])},GMAT[JJ,JJ]);

EMAT[1,1]:=diff(scm,u[j,n+1]);
EMAT[1,1+1]:=diff(scm,u[j+1,n+1]);
EMAT[JJ,JJ-1]:=diff(scm,u[j-1,n+1]);
EMAT[JJ,JJ]:=diff(scm,u[j,n+1]);

#print(GMAT);
#Given the form of GMAT and EMAT - evolve!!!

for ij from 2 to NN do
 print(`evolution`,ij);	

for i from 2 to JJ-1 do
MMAT[i,i]:=subs({u[j,n]=VEC[i],u[j-1,n]=VEC[i-1],
                 u[j+1,n]=VEC[i+1]},GMAT[i,i]);
MMAT[i,i+1]:=subs({u[j,n]=VEC[i],u[j+1,n]=VEC[i+1],
                   u[j-1,n]=VEC[i-1]},GMAT[i,i+1]);#print(4,");
MMAT[i,i-1]:=subs({u[j,n]=VEC[i],u[j-1,n]=VEC[i-1],
                   u[j+1,n]=VEC[i+1]},GMAT[i,i-1]);

EMAT[i,i]:=diff(scm,u[j,n+1]);
EMAT[i,i+1]:=diff(scm,u[j+1,n+1]);
EMAT[i,i-1]:=diff(scm,u[j-1,n+1]);
od;
#some extras!!!!!!!!!!!!!!!!!!!!
MMAT[1,1]:=subs({u[j,n]=VEC[1],u[j-1,n]=0,u[j+1,n]=VEC[2]},GMAT[1,1]);
MMAT[1,2]:=subs({u[j,n]=VEC[1],u[j-1,n]=0,u[j+1,n]=VEC[2]},GMAT[1,2]);
MMAT[JJ,JJ-1]:=subs({u[j,n]=VEC[JJ],u[j-1,n]=VEC[JJ-1],
                     u[j+1,n]=0},GMAT[JJ,JJ-1]);
MMAT[JJ,JJ]:=subs({u[j,n]=VEC[JJ],u[j-1,n]=VEC[JJ-1],u[j+1,n]=0},GMAT[JJ,JJ]);

EMAT[1,1]:=diff(scm,u[j,n+1]);
EMAT[1,1+1]:=diff(scm,u[j+1,n+1]);
EMAT[JJ,JJ-1]:=diff(scm,u[j-1,n+1]);
EMAT[JJ,JJ]:=diff(scm,u[j,n+1]);

#print(MMAT);print(EMAT);print(VEC);


#make the `product vector'
    VEC:=multiply(MMAT,VEC);
#evolve VEC
    VEC:=linsolve(EMAT,VEC);
#copy vector into SMAT
    for ii from 1 to JJ do SMAT[ij,ii]:=VEC[ii]; od;
  
od;


#matrixplot(SMAT,axes=BOXED);
#,labels=[`time`,`x`,`amplitude`]);

print(`mesh complete`);

end:

#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
#>


CROSS1:=proc(NN,JJ,MAT,pltvec,Nplt)
           local p,i,j,k,Q;

NP:=iquo(NN,Nplt);#print(NP);

k:=0;
#jumping thro time loop
for i from 1 by NP to NN do
k:=k+1;

L:=[];
#scanning thro space loop
for j from 1 to JJ do
l:=[[j,MAT[i,j]]];
L:=[op(L),op(l)];
od;

pltvec[k]:=plot(L,title=`waves at time intervals`);

od;

print(`plot structures in vector`);

end:

#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

with(plots):with(linalg):