USE stokes1

TYPEOF(zeta) zeta2old
IF PlotOn THEN GRAPHICSMODE
nexttime=time(NULL)+2
FILE risultati=CREATE("risultati")
DO psi(i)=y(i)∗EXP(-y(i)∗y(i)/tin); zstar(i)=psi(i) FOR ALL i
alfa=0.07; LOOP WHILE alfa<1.2-0.0000001
alfa=alfa∗1.1
sigma=(0.1+0.7∗I)∗alfa∗Re
t=tin
COMPLEX ampl=1
BaseFlow(); BuildMats()
T=A+sigma∗B; LUdecomp T
zeta=B∗psi
DO psi=T\zeta; zeta=B∗psi; zstar=zstar∗B; zstar=zstar/T FOR 5 TIMES
zetaold=zeta
dt=(tfin-tin)/nt
LOOP WHILE t<tfin-0.000001
  t=t+dt
  zeta2old=zetaold; zetaold=zeta
  eigenval
  COMPLEX c1=0
  DO c1=c1+zstar(i)∗[1.5∗zeta(i)-2∗zetaold(i)+0.5∗zeta2old(i)]/dt FOR i=zeta.LO TO zeta.HI
!  DO c1=c1+zstar(i)∗[zeta(i)-zetaold(i)]/dt FOR i=zeta.LO TO zeta.HI
  ampl=ampl∗EXP((sigma-c1)∗dt)
  IF t>tin+2.5∗dt THEN risultati Re∗sqrt(t), alfa∗sqrt(t), REAL(sigma/(I∗alfa∗Re)), -IMAG(sigma/(I∗alfa∗Re)), REAL(sigma/(I∗alfa∗Re)-c1/(I∗alfa∗Re)), -IMAG(sigma/(I∗alfa∗Re)-c1/(I∗alfa∗Re)), ampl
IF PlotOn THEN
  IF time(NULL)>=nexttime THEN
    REAL psimax=ABS[psi(1)]
    DO psimax=MAX(psimax,ABS[psi(i)]) FOR i=2 TO ny-1
    CLEARSCREEN; Plot(psi,8,psimax)
    WRITE "time=" t, "alfa=" alfa
    WRITE "ph.speed=" sigma/(I∗alfa∗Re) " (" (sigma-c1)/(I∗alfa∗Re) ")"
    WRITE "amplitude=" ampl
    nexttime=time(NULL)+2
  END IF
END IF
REPEAT LOOP
REPEAT LOOP

!(
FILE WriteProfile=CREATE("profile")
WriteProfile 0,0,0
DO WriteProfile y(i),ABS(I∗alfa∗psi(i)),ARG(I∗alfa∗psi(i)) FOR i=1 TO ny-1
WriteProfile
WriteProfile 0,0,0
DO WriteProfile (y(i)+y(i-1))/2, ABS([psi(i)-psi(i-1)]/[y(i)-y(i-1)]), ARG([psi(i)-psi(i-1)]/[y(i)-y(i-1)]) FOR i=2 TO ny-1
WriteProfile
WriteProfile 0,0,0
DO WriteProfile y(i),ABS(zstar(i)/(y(i+1)-y(i))∗2),ARG(zstar(i)/(y(i+1)-y(i))∗2) FOR i=1 TO ny-1
!)
IF PlotOn THEN TEXTMODE