	PROGRAM runge2

!********************************************************
!	Solves the ordinary differential equation dx(t)/dt=
!       -(sin(t³)+3t³ cos(t³))x(t)
!	with the heun's method and runge-kutta 4
!       Plot in Dislin
!       Compile with dlink -c runge2
!       Date: 2010
!********************************************************
	PARAMETER(T=3., DT=0.01, IDIM=T/DT)	
	DIMENSION TRAY(IDIM),XHEUN(IDIM),XRUNGE(IDIM)
	REAL FUNCX,TIME,XH
	Write(6,*) 'Gib x(t=0)'

	PI=2.*ASIN(1.)
        
!       Initial value of X
	Write(6,*) 'Initial value of X'
	READ(5,*) X

	XH=X

!       Calculating the K-Values for Runge-Kutta
	DO I=1, IDIM
	 TRAY(I)=DT*I	 
	 TIME=DT*I
	 XRUNGE(I)=X
	 XHEUN(I)=XH

	 XKEINS=FUNCX(X,TIME)*DT
	
	 XKZWEI=FUNCX(X+0.5*XKEINS,TIME)*DT
	
	 XKDREI=FUNCX(X+0.5*XKZWEI,TIME)*DT
	
         XKVIER=FUNCX(X+XKDREI,TIME)*DT
	
	 
    	 X=X+(XKEINS+2.*XKZWEI+2.*XKDREI+XKVIER)/6.
	 XH=XH+(XKEINS+2.*XKZWEI)/2.
C	 WRITE(6,*) TRAY(I)

	ENDDO
	
	CALL PLOT(TRAY,XHEUN,XRUNGE,IDIM,T)
	
		
	
	STOP
	END

C****************************************

	REAL FUNCTION FUNCX(X,TRAY,TIME)
	
	FUNCX=-(SIN(TIME**3)+3.*TIME**3*COS(TIME**3))*X
	RETURN 
	END 

C****************************************
!	X-t Diagramm

       SUBROUTINE PLOT(TRAY,XHEUN,XRUNGE,IDIM,T)
		
       DIMENSION TRAY(IDIM),XHEUN(IDIM),XRUNGE(IDIM),EXACT(IDIM)
       CHARACTER*30 CA
       CALL METAFL('XWIN')
       CALL DISINI
       CALL PAGERA
       CALL COMPLX
 
       CALL AXSPOS(450,1800)
       CALL AXSLEN(2200,1200)

       CALL NAME('X-axis','X')
       CALL NAME('Y-axis','Y')


       CALL GRAF(0.,T,0.,T/10.,-5.,20.,0.,5.)
     

       CALL COLOR('FORE') 
       CALL HEIGHT(50)
       CALL TITLE
       
!      WRITING EXACT SOLUTION TO EXACT(IDIM)
       DO I=1, IDIM
	EXACT(I)=EXP(-TRAY(I)*SIN(TRAY(I)**3))
       ENDDO
C       CALL INCMRK(-20)
C       CALL MARKER(21)
C       CALL HSYMBL(1)
       CALL COLOR('RED')
       CALL CURVE(TRAY,XHEUN,IDIM)
       CALL COLOR('GREEN')
       CALL CURVE(TRAY,XRUNGE,IDIM)
       CALL COLOR('BLUE')
       CALL CURVE(TRAY,EXACT,IDIM)
       
       CALL COLOR('FORE')
       CALL LEGINI(CA,3,10)
       CALL LEGTIT(' ')
       CALL LEGLIN(CA,'HEUN',1)
       CALL LEGLIN(CA,'R-K 4',2)
       CALL LEGLIN(CA,'EXACT',3)
       CALL LEGEND(CA,8)
       Write(6,*) 'globaler Diskretisierungsfehler Heun'
     & ,ABS(XHEUN(IDIM)-EXACT(IDIM))
       Write(6,*) 'globaler Diskretisierungsfehler Runge-Kutta 4'
     & ,ABS(XRUNGE(IDIM)-EXACT(IDIM))
       CALL DISFIN

   
       RETURN
       END

