Non-linear ODE Finite DIfference

    !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)
    !Matrix form:
    ![1     0            0     0               0]      [u_1]               [1]
    ![-1  1+2*delX*u_g2  0     0    ...        0]      [u_2]               [delX*u_g2^2]
    ![0    -1   1+2*delX*u_g3  0    ...        0]      [u_3]       =       [delX*u_g3^2]
    ![0     0           -1   1+2*delX*u_g4 ... 0]      [u_4]               [delX*u_g4^2]
    ![..........................................]      [...]               [...]
    ![0     0     0    ...   -1    1+2*delX*u_gN]_NxN  [u_N]_Nx1           [delX*u_gN^2]_Nx1
    ![A][U]=[B]
    !=> [U]=Inverse([A]).[B]
    !guess u_gi;i=1,2,...N and find [U] vector i.e. u_i;i=1,2,...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

    real A(50,50),U(50,50),B(50,50),U_guess(50,50)
    real AInverse(50,50)
    real error(50,50)
    real sumOfErrors,tolerance


integer N
integer i,j
real delX
N=4
delX=1/real(N-1)
U_guess=1 !fill all elements of initial guess vector with 1
tolerance=0.00001

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

10 A(1,1)=1
do i=2,N
    do j=1,N
        if(i.EQ.j) then
            A(i,j)=1+2*delX*U_guess(i,1)
        else if(i.EQ.j+1) then
            A(i,j)=-1
        else
            A(i,j)=0
        endif
    end do
    end do

    !populate the B vector
    B(1,1)=1
    do i=2,N
        B(i,1)=delX*U_guess(i,1)**2
    end do





    call Inverse(A,AInverse,N)
    call matrixMul(AInverse,B,U,N,N,N,1)

    !populate the error vector
    sumOfErrors=0
    do i=1,N
        error(i,1) =  abs(U(i,1)  -  U_guess(i,1))
        sumOfErrors=sumOfErrors+error(i,1)
    end do

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

    U_guess=U !update guess vector with solution of current iteration

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

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

    end program




subroutine matrixMul(A,B,C,arow,acol,brow,bcol)
integer arow,acol,brow,bcol
!To see how this works, try writing the elements of the product matrix in sigma notation
dimension A(50,50),B(50,50),C(50,50) !CAUTION: This declaration is crucial
do k=1,arow !k runs from 1 to arow
    do l=1,bcol !and l from 1 to bcol because that's the output matrix dimension and we _have_ to populate each element of it
        sum=0
        do i=1,acol !this loop is the essence of the aforementioned sigma notation for the (k,l) element of the product matrix
            sum=sum+A(k,i)*B(i,l)
        end do
        C(k,l)=sum
    end do
end do
return
end subroutine matrixMul

subroutine Inverse(X,XI,size)
dimension X(50,50),XI(50,50),originalX(50,50)
integer rowindex,colindex,size,tmpi,tmpj,i,j
real pivotelement

originalX=X

do i=1,size
    do j=1,size
        if(i.EQ.j) then
            XI(i,j)=1
        else
            XI(i,j)=0 !this is for good measure coz by default, XI may have been changed by the caller itself though it is 0 by default
        endif
    end do
end do


do rowindex=1,size !loop over rows of X

    pivotelement=X(rowindex,rowindex) !this is crucial!

    !divide the row by pivot element
    do colindex=1,size
        X(rowindex,colindex)=X(rowindex,colindex)/pivotelement
    end do
    do colindex=1,size
        XI(rowindex,colindex)=XI(rowindex,colindex)/pivotelement
    end do


    !perform row operations to make non-pivot elements of this column zero
    do i=1,size !i gives the row number being manipulated
        pivotelement=X(i,rowindex) !pivot element for the currently manipulated row. this i crucial

        if(i.NE.rowindex) then !row operations are for non-current rows only

            do j=1,size !column index for ith row of X matrix
                X(i,j)=X(i,j)-pivotelement*X(rowindex,j)
            end do
            do j=1,size !column index for ith row of XI matrix
                XI(i,j)=XI(i,j)-pivotelement*XI(rowindex,j)
            end do


        endif
    end do



end do

!restore X to original X since upon inversion completion, X has been reduced to an Identity matrix. good practice
X=originalX

end subroutine