Initializing and Solving matrices

here is my whole program. It will be used to solve bvp's (boundary Value Problems) using the finite difference method. The problem i have is initializing the matrix that needs to be solved and then having it printed out. I've spent a lot of time on it today and i hope to figure it out later on. If anyone has any suggestions on how to improve the code so that i can actually get the matrix initialized and the solution printed out, all help will be greatly appreciated. By the way, the formulas used are correct
and the only thing that is having problems is initializing the A matrix and printing out the X matrix, which is the solution.

Thanks

Code:

---

//All Header Files Used
#include <iostream>
#include <cstdlib>
#include <string>
#include <cmath>
#include <fstream>
#include <cctype>
#include <climits>

//initialising all matrices and vectors to the size of 100
double A[100][100],L[100][100],f[100], d[100], X[100],U[100][100];

//declaring functions for use in the program
void Finite_Difference(int n,double a,double b,double c,double dx,double A[ ][100]);
void LU(int n, double A[][100],double L[][100]);
void f_sub(int n,double L[ ][100],double f[100],double g[100]);
void b_sub(int n,double U[ ][100],double g[100],double X[100]);

using namespace std;

//Function Definition to Calculate the Finite Difference Matrix
void Finite_Difference(int n,double a,double b,double c,double dx, double A[ ][100])
{
       
        for(int i=0;i<n;i++)
        {
                for(int j=0;j<n;j++)
                {
                        if (i == j)
                        {
                                A[i][j]=((2*a) - c*(pow(dx,2)));
                        }

                        else if( j == (i+1))
                        {
                                A[i][j] = (-a + (b/2)*dx);
                        }

                        else if (j == (i-1))
                        {
                                A[i][j] = (-a - (b/2)*dx);
                        }

                        else
                                A[i][j] = 0;
                }
        }
}

//Function Definition to Factor the A matrix in to the [L][U] matrices
void LU(int n,double U[ ][100],double L[ ][100] )
{
        for(int j=0;j<(n-1);j++)
        {
                for(int i=j+1;i<n;i++)
                {
                        L[i][j]= (U[i][j]/U[j][j]);

                        for(int k=j+1;k<n;k++)
                        {
                                U[i][k]=U[i][k]-(L[i][j]*U[j][k]);
                        }
                }
        }

        for(int i=0;i<n;i++)
        {
                for(int j=0;j<n;j++)
                {
                        if(i==j)
                                L[i][j]=1;
                }
        }
}

//Function Definition for Forward Substitution
void f_sub(int n,double L[ ][100],double f[100],double g[100] )
{
        g[0]=f[0];
        for(int i=1;i<n;i++)
        {
                g[i]=f[i];
                for(int k=0;k<i;k++)
                {
                        g[i]=g[i]-(L[i][k]*g[k] );
                }
        }
}


///Function Definition for Backward Substitution
void b_sub(int n,double U[ ][100],double g[100],double X[100] )
{

        X[n-1]= g[n-1]/U[n-1][n-1];

        for(int i=n-2;i>=0;i--)
        {
                X[i]=g[i];
                for(int k=i+1;k<n;k++)
                {
                        X[i]=X[i]-(U[i][k]*X[k]);
                }

                X[i]=X[i]/(U[i][i] );
        }
}


int main()
{

        //initialises a boolean so the user can choose to quit or continue another trial
        bool again = true;

        while(again)
        {
                cout<<" Welcome to the Boundary Value Problem Solver"<<endl;
                cout<<" This Solver Will be Using the Finite Difference Method"<<endl;
                cout<<"        The Differential Equation to be Solved is in the Form of:"<<endl;
                cout<<endl;
                cout<<"                        ay'' + b' + cy + x = 0"<<endl;
                cout<<" With Boundary Conditions y(0) = d & y(20) = e"<<endl;
                double const1,const2,const3,bound1,bound2,dx;
                bool v1 = false;
                while(!v1)
                {
                        cout<<"Please Enter the Constant 'a':\n";
                        cin>>const1;

                        //Determined whether the input is valid or invalid
                        if(cin.fail())
                        {
                                cin.clear();
                                cin.ignore(INT_MAX, '\n');
                                cout << "The input was invalid. Please try again" << endl;
                        }
                       
                        else
                        {
                                v1 = true;
                        }
                }

                bool v2 = false;
                while(!v2)
                {
                        cout<<"Please Enter the Constant 'b':\n";
                        cin>>const2;

                        //Determined whether the input is valid or invalid
                        if(cin.fail())
                        {
                                cin.clear();
                                cin.ignore(INT_MAX, '\n');
                                cout << "The input was invalid. Please try again" << endl;
                        }
                       
                        else
                        {
                                v2 = true;
                        }
                }

                bool v3 = false;
                while(!v3)
                {
                        cout<<"Please Enter the Constant 'c':\n";
                        cin>>const3;

                        //Determined whether the input is valid or invalid
                        if(cin.fail())
                        {
                                cin.clear();
                                cin.ignore(INT_MAX, '\n');
                                cout << "The input was invalid. Please try again" << endl;
                        }
                       
                        else
                        {
                                v3 = true;
                        }
                }

                bool v4 = false;
                while(!v4)
                {
                        cout<<"Please Enter the Boundary Condtion at y(0):\n";
                        cin>>bound1;

                        //Determined whether the input is valid or invalid
                        if(cin.fail())
                        {
                                cin.clear();
                                cin.ignore(INT_MAX, '\n');
                                cout << "The input was invalid. Please try again" << endl;
                        }
                       
                        else
                        {
                                v4 = true;
                        }
                }

                bool v5 = false;
                while(!v5)
                {
                        cout<<"Please Enter the Boundary Condition at y(20):\n";
                        cin>>bound2;

                        //Determined whether the input is valid or invalid
                        if(cin.fail())
                        {
                                cin.clear();
                                cin.ignore(INT_MAX, '\n');
                                cout << "The input was invalid. Please try again" << endl;
                        }
                       
                        else
                        {
                                v5 = true;
                        }
                }

                bool v6 = false;
                while(!v6)
                {
                        cout<<"Please Enter the Interval Size:\n";
                        cin>>dx;

                        //Determined whether the input is valid or invalid
                        if(cin.fail())
                        {
                                cin.clear();
                                cin.ignore(INT_MAX, '\n');
                                cout << "The input was invalid. Please try again" << endl;
                        }
                        else if(dx<=0)
                        {
                                cout<< "The Interval Must be Greater than 0"<<endl;
                                cout<< "Please Try Again."<<endl;
                        }
                        else
                        {
                                v6 = true;
                        }
                }

                double dim = ((bound2-bound1)/dx) - 1;
                int size = (int)dim;

                //Initialises the A,L, and U matrices to 0 everywhere
                for(int i=0;i<size-1;i++)
                {
                        for(int j=0;j<size-1;j++)
                        {
                                A[i][j]=0;
                                L[i][j]=0;
                                U[i][j]=0;
                        }
                }

                //creates the Finite_Difference matrix
                Finite_Difference(size,const1,const2,const3,dx,A);

                //Create ouput text file
                ofstream output_data;
                output_data.open("output.txt");

                //Check to see if output file has opened correctly
                if(output_data.fail())
                {
                        cout<<"Output file opening failed.\n";
                        exit(1);
                }

                //creates the Right-Hand Vector
                int xi = 0;
                        for( i=0;i<size;i++)
                        {
                                xi+=dx;
                                if (i==0)
                                {
                                        f[i] = bound1 + (xi*(pow(dx,2)));
                                }

                                else if (i == (size -1))
                                {
                                        f[i] = bound2 +(xi*(pow(dx,2)));
                                }

                                else
                                {
                                        f[i] = xi*(pow(dx,2));
                                }
                        }
                       
                cout<<endl;

                LU(size,A, L);
                f_sub( size,L, f, d);
                b_sub( size,A, d, X);

                output_data<<"The Solution is:"<<endl;
                output_data<<"                x                                y"<<endl;
                output_data<<"        _________                __________"<<endl;
                output_data<<"  "<<0<<"                        "<<bound1<<endl;

                cout<<"The Solution to this Problem can be found in the file <output.txt>:\n";
                int ex = dx;
                for( i=0;i<size;i++)
                {
                        ex+=dx;
                        output_data<<"  "<<ex<<"                        "<<X[i]<<endl;
                }
                output_data<<"  "<<20<<"                        "<<bound2<<endl;

                cout<<"Would you like to Run This Program Again?(y/n)"<<endl;
                string encore;
                cin>>encore;

                //restarts another analysis
                if(encore == "y" || encore == "Y")
                {
                        //clears the screen
                        system("cls");
                        again = true;
                }

                //quits the program
                else if(encore == "n" || encore == "N")
                {
                        again = false;
                }

                //Check whether the input is valid
                else
                {
                        cout<<"You Have not Entered a Proper Argument. Quiting Anyways"<<endl;
                        again = false;
                }
        }


        return 0;
}

---

oh and the method of solving the matrix used in the program is by LU decomposition.