Adams-Bashforth 4 step method diverging

Hi - student trying to write a program to solve a trivial 1st order DE (dy/dx=-xy), the idea being that once it works it can be used for different forms of f(x,y). I use Euler to get the starting point, then A-B's 2 step formula for the next 2 starting points.
It is when I next apply the 4-step formula that I get obviously wrong values for Y(n+1) (using debug). I've looked at this every which way until my head hurts, can't see what I'm doing wrong. Help would be nice javascript:emoticon(':shock:')

program AB4step
  	!Summary: Apply Adams-Bashforth 4 step method to dy/dx = -xy given y(0)=1
    implicit none
    integer, parameter  :: dp = selected_real_kind(15, 307)
	real(kind=dp)	   	:: xN,yN,xNplus1,yNplus1,xNminus1,yNminus1,exact,diff
	real(kind=dp)	   	:: fN, fNminus1, fNminus2, fNminus3, FUNC
    real, parameter     :: lim=4.0
    real                :: h
    integer				:: iter, numsteps
	integer, parameter  :: n=60		
    character(len=n )   :: stars=REPEAT('*',n) 
   
    !----- header ----                                   
	print *, ' ' 
	print 10, 'Program AB-4step starts, integrating from 0.0 to ', lim
    print *, 'Enter the number of steps to use (zero or negative to stop)'
    read  (*,*) numSteps            ! ensure integer number os steps, adjust h accordingly
    print *, stars
 	print *, ' ' 
	print 20, 'Iter.','X(n+1)','F(n)    ', 'Y(n+1)      ', 'Exact   ','Diff    '
	!----- end header ----                                   

	do while (numSteps.GT.0)    	! allows trials of diff. values of step size h
	!----- Initialise ----                                   
        h=lim/numSteps
        yNminus1=1.0        ! y(0)              starting point, n=0
        xNminus1=0.0        ! x=0
        fNminus2=0.0
        fNminus1=FUNC(xNminus1, yNminus1)
        yN=yNminus1+(h*fNminus1)   ! using Eulers 1-step to get 2nd starting value, n=1
        xN=xNminus1+h

    !----- Calculate next 2 points using AB 2step (n = 2,3) --------    
        Do iter=1,2                         
            fN=FUNC(xN,yN)
            exact=exp(-(xN**2)/2.0)
            diff=exact-yN
            print 50, iter, xN, fN, yN, exact, diff

            yNplus1=yN+(h*((1.5*fN)-(0.5*fNminus1)))
            xNplus1=xN+h
            yNminus1=yN
            xNminus1=xN
            fNminus3=fNminus2
            fNminus2=fNminus1
            fNminus1=fN
            yN=yNplus1
            xN=xNplus1         
     	end do 
    !----- Calculate rest of points using AB 4 step (n = 4 onward --------    
        Do iter=iter,numsteps                          
            fN=FUNC(xN, yN)
            exact=exp(-(xN**2)/2.0)
            diff=exact-yN
            print 50, iter, xN, fN, yN, exact, diff

            yNplus1=yN+((h/24.0)*((55.0*fN)-(59.0*fNminus1)+(37.0*fNminus2)-(9.0*fNminus3)))
            xNplus1=xN+h
            yNminus1=yN
            xNminus1=xN
            fNminus3=fNminus2
            fNminus2=fNminus1
            fNminus1=fN
            yN=yNplus1
            xN=xNplus1         
     	end do 
        print *, 'The step size used above was ',h
        print *, 'Enter another number of steps value (zero or negative to stop)'
        read  (*,*) numSteps
    end do
    
! admin     
    print *, 'User selected end'                                  
  	print *, stars  
10  format (A, F8.3)
20  format (A, 5X, 5(A,8X))
50  format (I2, 5X, F8.4, 2X, 4(F15.9, 3X))
  	end  program AB4step

!---------Functions -----------
    function FUNC (x, y)
    integer, parameter  		:: dp = selected_real_kind(15, 307)
    real(kind=dp)	   			:: x, y, FUNC
        FUNC=-x*y
    end function

Typically students like you who do not care what they post do not get much attention on any forums

? What have I done wrong?

Quoted from ognik ? What have I done wrong?

Here are a few hints.

Your code is pasted in-line, causing the indentation to get lost, and making the code unpleasant to read. Use code.../code tags around your code fragments (surrounded by brackets).

Your code is longer than the forum limit, so the end part is trimmed (have you noticed this?), making the program incomplete and unusable. Keep them short, or place longer files elsewhere and provide a link.

You state that the result is wrong, but do not state what you obtained, what you think the correct result is, and how close you expect the two to be to each other. Explain your reasoning.

Here is a fragment of code for you to work with. Add declarations before it, complete the AB-4 part, add output statements after the fragment, and try. After you get the code to run and produce satisfactory results, you can consider using x, y and f values at only four x-levels instead of using arrays (these are convenient and make the code easy to read, but could be considered wasteful of memory).

!
! Initial condition
x(0)=0.; y(0)=1.; h=0.01
!
! Euler step
!
f(0)=fn(x(0),y(0))
x(1)=h; y(1)=y(0)+h*f(0)
f(1)=fn(x(1),y(1))
!
! AB-2 steps
!
DO i=2,3
   x(i)=x(i-1)+h
   y(i)=y(i-1)+(3.*f(i-1)-f(i-2))*0.5*h
   f(i)=fn(x(i),y(i))
END DO
!
! AB-4 steps
!
DO i=4,N

My apologies and thanks for the advice. I did review the code before submitting - and it was fine in review, however after submitting it lost the last few lines, I will remember the length limit.

Actually my program turned out to be correct - I was simply using too large a steps-size, which appears to cause A-B 4th to oscillate. But thanks for the advice - I am learning Fortan and clearly using arrays is better than my clumsy names .