program temp
c      Analysis of vanderpol oscillator. A*Sin(om*t) is the excitation,
c      mu*(x^2-1) is the damping coefficient.
       integer i,j,nsteps,niter,printstep
       real*8 a,om,mu,x0,v0,xn,xnp1k,xmk,vn,vm,vnp1,xmdelta,
     & tn,tnp1,xnp1delta,tdelta,fdelta1,fdelta2,tprint,
     & pi,m11,m12,m21,m22,detm,tend,uf(300000),vf(300000),
     & time(300000),energy(300000)
       open(11,file='shape',status='unknown')
c      ************************************************
       pi=4.0d0*datan(1.0d0)
       a=0.0d0
       om=3.0d0
       mu=1.0d1
       x0=1.0d0
       v0=0.0d0
       tend=100.0d0
       niter=10
       tdelta=0.001d0
       tprint=0.01d0
c      prints after every tprint times
       xmdelta=0.0d0
       xnp1delta=0.0d0
       nsteps=nint(tend/tdelta)
       printstep=nint(tprint/tdelta)
       if(printstep.lt.1)printstep=1
       if(nint(1.0d0*nsteps/printstep).gt.300000)then
        write(11,*) 'Increase size of arrays uf, vf etc.'
        stop
       endif
       print*,'nsteps,niter ',nsteps,niter
       time(1)=0.0d0
       xn=x0
       vn=v0
       uf(1)=x0
       vf(1)=v0
       energy(1)=(v0**2+x0**2)/2.0d0
       xmk=x0
       xnp1k=x0
       do 1 i=2,nint(1.0d0*nsteps/printstep)+1
        time(i)=printstep*(i-1)*tdelta
1      continue
       do 100 i=1,nsteps
       tn=(i-1)*tdelta
       tnp1=i*tdelta
       do 10 j=1,niter

      m11= tdelta**2*(140 + 8*mu*vn*(-8*xmk + xn)) + 8*mu*tdelta*
     & (72*xmk**2 - 7*(5 + xn**2) + 2*xnp1k**2 - 4*xmk*(7*xn + xnp1k))

      m12=420 + 2*mu*tdelta**2*vn*(-xn + xnp1k) + 2*mu*tdelta*(-35 
     & - 2*xn*xnp1k + 21*xnp1k**2 + 8*(-xmk**2 + xn**2 + 2*xmk*xnp1k))

      m21=-3360 + tdelta**2*(140 + 8*mu*vn*(-8*xmk + xnp1k)) + 
     & 8*mu*tdelta*(35 + 24*xmk**2 - 3*xn**2 + 6*xn*xnp1k 
     & - 14*xnp1k**2 - 4*xmk*(7*xn + 5*xnp1k)) 

       m22=1260 + 2*mu*tdelta*(-175 - 40*xmk**2 - xn**2 + 
     & 8*xmk*(3*xn - 14*xnp1k) + 76*xn*xnp1k + 228*xnp1k**2) + 
     & tdelta**2*(70 + 2*mu*vn*(4*xmk - xn + 25*xnp1k))

        detm=m11*m22-m12*m21

      fdelta1=-(420*(-xn + xnp1k) + tdelta**2*(140*xmk + 70*xn + 
     & mu*(-32*vn*xmk**2 + 8*vn*xmk*xn + 25*vn*xn**2 - 2*vn*xn*xnp1k 
     &+ vn*xnp1k**2)) + tdelta*(-420*vn + 2*mu*(96*xmk**3 - 26*xn**3 
     & - 28*xmk*(5 + xn**2) + xn*(175 - xnp1k**2) + 7*xnp1k*
     & (-5 + xnp1k**2) + 8*(xn**2*xnp1k + xmk*xnp1k**2 - 
     & xmk**2*(7*xn + xnp1k)))) + (420*A*(-(om*tdelta*Cos(om*tn)) 
     & - Sin(om*tn) + Sin(om*tnp1)))/om**2)

        fdelta2=-(420*(-8*xmk + 5*xn + 3*xnp1k) + tdelta**2*(140*xmk 
     & + 70*xnp1k + mu*(-32*vn*xmk**2 + vn*xn**2 + 8*vn*xmk*xnp1k - 
     & 2*vn*xn*xnp1k + 25*vn*xnp1k**2)) + tdelta*(420*vn + 2*mu*
     & (32*xmk**3 + 35*xn - 5*xn**3 - (175 + xn**2)*xnp1k + 
     & 38*xn*xnp1k**2 + 76*xnp1k**3 - 8*xmk**2*(7*xn + 5*xnp1k) + 
     & 4*xmk*(35 - 3*xn**2 + 6*xn*xnp1k - 14*xnp1k**2))) + 
     & (420*A*(om*tdelta*Cos(om*tnp1)+Sin(om*tn)-Sin(om*tnp1)))/om**2)

       xmdelta=(m22*fdelta1-m12*fdelta2)/detm
       xnp1delta=(-m21*fdelta1+m11*fdelta2)/detm
       xmk=xmk+xmdelta
       xnp1k=xnp1k+xnp1delta
       if(dsqrt(xmdelta**2+xnp1delta**2).lt.1.0d-10)go to 20
10     continue
       print*,'Failed to converge'
       stop
20     continue
       vm=(4*xmk + xnp1k -5*xn - tdelta*vn)/(2*tdelta)
       vnp1=4*(xnp1k+xn-2*xmk)/tdelta + vn
       vn=vnp1

       if(i.ge.printstep)then
       if(mod(i,printstep).eq.0)then
       uf(nint(1.0d0*i/printstep)+1)=xnp1k
       vf(nint(1.0d0*i/printstep)+1)=vnp1
       energy(nint(1.0d0*i/printstep)+1)=(vnp1**2+xnp1k**2)/2.0d0
       endif
       endif
       xn=xnp1k
100    continue

       do 200 i=1,nint(1.0d0*nsteps/printstep)+1
       write(11,500) time(i),uf(i),vf(i),energy(i)
200    continue
500    format(1x,e14.7,1x,e14.7,1x,e14.7,1x,e14.7)
       end