c
c     Numerical Analysis:
c     The Mathematics of Scientific Computing
c     D.R. Kincaid & E.W. Cheney
c     Brooks/Cole Publ., 1990
c
c     Section 8.6
c
c     Taylor series method (order 3) for solving a
c     system of ordinary differential equations
c
c
c     file: taysys.f
c
      data  t/1.0/, x/3.0/, y/1.0/, h/-0.1/, m/30/   
c
      print *
      print *,' Taylor series method (order 3) for a system'
      print *,' Section 8.6, Kincaid-Cheney'
      print *
      print 3,'k','t','x','y' 
      print 4,0,t,x,y
c
      do 2 k = 1,m
        x1 = x + y**2.0 - t**3.0      
        y1 =  y + x**3.0 + cos(t)   
        x2 = x1 + 2.0*y*y1 - 3.0*(t**2.0)    
        y2 =  y1 + 3.0*x1*(x**2.0) - sin(t) 
        x3 = x2 + 2.0*y*y2 + 2.0*(y1**2.0) - 6.0*t    
        y3 = y2 + 6.0*x*(x1**2.0) + 3.0*(x**2.0)*x2 -cos(t)   
        x = x + h*(x1 + (h/2.0)*(x2 + (h/3.0)*(x3)))
        y = y + h*(y1 + (h/2.0)*(y2 + (h/3.0)*(y3)))
        t = t + h       
        print 4,k,t,x,y       
 2    continue    
c
 3    format(a6,a9,2a15)
 4    format(1x,i5,2x,3(e13.6,2x))
      stop
      end