Non-linear ODE FiniteDIfference algebraic iterative solution

    !Discretization of du/dx + u^2 = 0 for 0<=x<=1, discretized to N-1 parts(so that x and u range from x_1 to x_N)
    !=> (u_i - u_(i-1))/delX + u_i^2 = 0
    !let u_i = u_gi + delu_i where, u_gi = guess for u_i and delu_i = error in guess
    !so, u_i^2 = u_gi^2 + 2*delu_i*u_gi + delu_i^2. Neglecting delu_i^2 for being very small,
    !u_i^2 = 2*u_gi*delu_i + u_gi^2 = 2*u_gi*(u_i - u_gi) + u_gi^2 = 2*u_i*u_gi - u_gi^2
    !So, linearized(in u_i) finite difference equation:
    !(2*delX*u_gi + 1)*u_i - u_(i-1) - delX*u_gi^2 = 0 ; i=2,3,4,..,N
    !and u_1 = 1 (Boundary Condition)
    !guess u_gi;i=1,2,...N and find u_i;i=1,2,...N as following:
    !u_i = (u_(i-1) + delX*u_gi^2)/(1 + 2*delX*u_gi) ; i=2,3,...N
    !update u_gi with the recently computed u_i and continue iterating till Sum(|u_i - u_gi|);i=1,2,...N is less than tolerance

    implicit none

    integer N,i,j
    real delX
    real U(50),U_guess(50),error(50)
    real sumofErrors,tolerance

    tolerance=0.00001
    N=4
    delX=1/real(N-1)

    write(*,*) "N = ",N

    U_guess=1 !populate initial guess

10  U(1)=1 !boundary condition: u_1=1, never changes
    do i=2,N
        U(i)=( U(i-1) + delX*U_guess(i)**2 ) / ( 1  +  2*delX*U_guess(i) )
    end do

    error=abs(U-U_guess) !populate error vector

    sumofErrors=sum(error)

    write(*,*) "U guess = ",(U_guess(i),i=1,N)
    write(*,*) "U       = ",(U(i),i=1,N)
    write(*,*)

    U_guess=U !update guess vector with current output

    if(sumofErrors.GT.tolerance) then
        goto 10
    endif

    open(5,file="output.txt")
    do i=1,N
        write(5,*) U(i)
    end do
    close(5)

end program