USE stokes1

ASK alfa
IF PlotOn THEN GRAPHICSMODE

!(
DO zeta(i)=EXP(-y(i)) FOR ALL i 
!)
zeta=0

zetaold=zeta

TYPEOF(zeta) zeta2old
t=tin; INTEGER it=0
nexttime=time(NULL)+1
FILE energy=CREATE("march-10")

BOOLEAN nonancora=TRUE

LOOP WHILE it<nt
  it=it+1
  dt=tin+(tfin-tin)∗it∗(0.5∗it+nt0)/nt/(0.5∗nt+nt0)-t
  t=t+dt
  zeta2old=zetaold; zetaold=zeta
  BaseFlow(); BuildMats()
  T=A+1.5/dt∗B
  DO zeta(i)=(2∗zetaold(i)-0.5∗zeta2old(i))/dt FOR ALL i

  IF (t>0.1) AND nonancora THEN
  zeta(1) = zeta(1) - T(1,-1)/dt
  zeta(2) = zeta(2) - T(2,-2)/dt 
  END IF

  LUdecomp T

  psi=T\zeta; zeta=B∗psi
  IF (t>0.1) AND nonancora THEN
  nonancora=FALSE
  zeta(1) = zeta(1) + B(1,-1)/dt 
  zeta(2) = zeta(2) + B(2,-2)/dt 
  END IF
  sigma=SUM psi(i)|[1.5∗zeta(i)-2∗zetaold(i)+0.5∗zeta2old(i)]∗[y(i+1)-y(i-1)]/2 FOR ALL i
  e=Energy(psi)
  IF(e#0)THEN sigma=sigma/e/dt
IF PlotOn THEN
  IF time(NULL)>=nexttime THEN
    REAL psimax=1.E-10
    DO psimax=MAX(psimax,ABS(psi(i))) FOR i=1 TO ny-1
    CLEARSCREEN
    WRITE "time="t, "Energy="Energy(psi), "sigma="sigma
    WRITE "ph.speed=" sigma/(I∗alfa∗Re)
    Plot(psi,8,psimax)
    nexttime=time(NULL)+1
  END IF
  WRITE TO energy Re∗sqrt(t), alfa∗sqrt(t), REAL(sigma/(I∗alfa∗Re)), -IMAG(sigma/(I∗alfa∗Re)) 
END IF
REPEAT LOOP

FILE Profile=CREATE("profile")
WRITE TO Profile 0,0,0
DO WRITE TO Profile y(i),ABS(I∗alfa∗psi(i)),ARG(I∗alfa∗psi(i)) FOR i=1 TO ny-1
WRITE TO Profile
WRITE TO Profile 0,0,0
DO WRITE TO Profile (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