C     Euler method for numerical solution for differential equations 
      dimension x(0:100),y(0:100),v(0:100)
      n=100
      call difeq(x,y,v,n)
      call init
      call plotd1(x,y,n)
      call plotd2(y,v,n)
      call fin
      stop
      end
c
      subroutine difeq(x,y,v,n)
      dimension x(0:100),y(0:100),v(0:100)
       a=1.0
       om=1.0
       x(0)=0.d0
       y(0)=1.0d0
       v(0)=a*om*sin(om*x(0))
       dx=2.0*3.14/dble(n)
       do 10 i=1,n
        x(i)=dble(i)*dx
        v(i)=a*om*sin(om*x(i))
        y(i)=y(i-1)-dx*v(i)
c        write(*,*) x(i),y(i)
   10  continue
      return
      end
c
      subroutine plotd1(x,y,n)
      dimension x(0:100),y(0:100),x1(0:100),y1(0:100)
      call xyworld(-0.2, -1.2, 2.2, 1.2)
      call viewport(0.05, 0.05, 0.45, 0.45)
      call frame1
      xmax=6.28
      ymax=1.0
      do 20 i=0,n
        x1(i)=x(i)/xmax*2.0
        y1(i)=y(i)/ymax
   20  continue
      call linety(1)
      call plot(x1(0),y1(0),3)
      do 30 i=1,n
        call plot(x1(i),y1(i),2)
   30 continue
      call stroke
      do 41 i=0,n
       call square(x1(i),y1(i),0.01,0.0,1.0,1)
   41 continue
      return
      end
      subroutine plotd2(y,v,n)
      dimension y(0:100),v(0:100),x1(0:100),y1(0:100)
      call init
      call xyworld(-1.2, -1.2, 1.2, 1.2)
      call viewport(0.55, 0.05, 0.95, 0.45)
      call frame2
      ymax=1.0
      vmax=1.0
      do 20 i=0,n
        x1(i)=y(i)/ymax
        y1(i)=v(i)/vmax
   20  continue
      call linety(1)
      call plot(x1(0),y1(0),3)
      do 30 i=1,n
        call plot(x1(i),y1(i),2)
   30 continue
      call stroke
      do 41 i=0,n
       call square(x1(i),y1(i),0.01,0.0,1.0,1)
   41 continue
      call stroke
      return
      end
      subroutine frame1()
       call arrowb( -0.2,  0.0, 2.4, 0.0, 0.02, 0.0, 1.0)
       call arrowb(  0.0, -1.1, 0.0, 1.2, 0.02, 0.0, 1.0)
       call textx(2.4,-0.02,1,'x')
       call textx(2.0,0.0,4,'6.28')
       call textx(1.0,0.0,4,'3.14')
       call texty(-0.02,1.2,1,'y')
       call texty(0.0, 1.0, 1,'1')
       call texty(0.01, 0.5, 3,'0.5')
       call texty(0.0,-1.0, 2,'-1')
       call texty(0.01,-0.5, 4,'-0.5')
       call text(-0.2,-0.2,1,'O')
       n=10
       do 10 i=1, n
        x=0.2*i
        if(i/5*5.eq.i)  then 
          dy=0.02
        else 
          dy=0.01
        end if
        call line(  x, -dy,  x, dy, 0.0, 1.0, 1)
   10  continue
       do 20 i=1, n
        y=0.1*i
        if(i/5*5.eq.i)  then 
          dx=0.02
        else 
          dx=0.01
        end if
        call line( -dx,  y, dx,  y, 0.0, 1.0, 1)
        call line( -dx, -y, dx, -y, 0.0, 1.0, 1)
   20  continue
      return
      end
      subroutine frame2()
       call arrowb( -1.1,  0.0, 1.2, 0.0, 0.02, 0.0, 1.0)
       call arrowb(  0.0, -1.1, 0.0, 1.2, 0.02, 0.0, 1.0)
       call textx(1.2,-0.02,1,'x')
       call textx(1.0,0.0,1,'1')
       call textx(0.5,0.0,3,'0.5')
       call textx(-1.0,0.0,2,'-1')
       call textx(-0.5,0.0,4,'-0.5')
       call texty(-0.02,1.2,1,'y')
       call texty(0.0, 1.0, 1,'1')
       call texty(0.01, 0.5, 3,'0.5')
       call texty(0.0,-1.0, 2,'-1')
       call texty(0.01,-0.5, 4,'-0.5')
       call text(-0.2,-0.2,1,'O')
       n=10
       do 10 i=1, n
        x=0.1*i
        if(i/5*5.eq.i)  then 
          dy=0.02
        else 
          dy=0.01
        end if
        call line(  x, -dy,  x, dy, 0.0, 1.0, 1)
        call line( -x, -dy, -x, dy, 0.0, 1.0, 1)
   10  continue
       do 20 i=1, n
        y=0.1*i
        if(i/5*5.eq.i)  then 
          dx=0.02
        else 
          dx=0.01
        end if
        call line( -dx,  y, dx,  y, 0.0, 1.0, 1)
        call line( -dx, -y, dx, -y, 0.0, 1.0, 1)
   20  continue
      return
      end