Beregning virker

#include <cmath>
#include <cstdio>
#include <fstream>
#include <iostream>
#include "nr3.h"
#include "svd.h"
#include "utilities.h"

using namespace std;
void message(string input) {
    auto sbar = "----------------------------------\n";
    cout << sbar << input << '\n' << sbar;
}
double F(double x, double y, double d) {
    return (0.5 / pow(d * d + (x - y) * (x - y), 3.0 / 2.0));
}
double I1(VecDoub points, double x, double alimit, double h, double d) {
    int N = points.size() / 2;
    double res = 0;
    double y = alimit;
    for (int i = 0; i < N; i++, y = y + h) {
        if (i == 0 || i == N - 1) {
            res += 0.5 * F(x, y, d) * points[i];
        } else {
            res += F(x, y, d) * points[i];
        }
        cout << h << endl;
    }
    return res;
}

int main() {
    // JOHNS SETION
    int N = 4;
    // Problem size N
    {
        MatDoub A(2 * N + 2, 2 * N + 2);
        VecDoub b(2 * N + 2);
        const double eps1 = 0.8;
        const double eps2 = 0.6;
        const double T1 = 1000.0;
        const double T2 = 500.0;
        const double d = 1.0;
        const double w = 1.0;
        const double boltz = 1.712 * pow(10, -9);
        double blimit = 0.5 * w;

        double alimit = -0.5 * w;
        double h = (blimit - alimit) / N;
        double k1 = (1 - eps1);
        double k2 = (1 - eps2);
        for (int i = 0; i < N + 1; i++) {
            A[i][i] = -1;
            b[i] = -eps1 * boltz * pow(T1, 4);
        }
        for (int i = N + 1; i < 2 * N + 2; i++) {
            A[i][i] = -1;
            b[i] = -eps2 * boltz * pow(T2, 4);
        }
        {
            double y = 0;
            double x;
            int yy = 0;
            int xx = 0;
            for (yy = 0, y = alimit; yy < N + 1; yy++, y = y + h) {
                for (xx = N + 1, x = alimit; xx < 2 * N + 2; xx++, x = x + h) {
                    // Insert correct value at edges
                    if (yy == 0 || yy == N) {
                        A[xx][yy] = ((h * k2) / 2) * F(x, y, d);
                    } else {
                        A[xx][yy] = (h * k2) * F(x, y, d);
                    }
                }
            }
        }

        {
            double y = 0;
            double x;
            int yy = 0;
            int xx = 0;
            for (xx = 0, y = alimit; xx < N + 1; xx++, y = y + h) {
                for (yy = N + 1, x = alimit; yy < 2 * N + 2; yy++, x = x + h) {
                    // Insert correct value at edges
                    if (yy == N + 1 || yy == 2 * N + 1) {
                        A[xx][yy] = ((h * k1) / 2) * F(x, y, d);
                    } else {
                        A[xx][yy] = (h * k1) * F(x, y, d);
                    }
                }
            }
        }
        VecDoub xVec(2 * N + 2);

        SVD svddcmp(A);
        svddcmp.solve(b, xVec);
        // Calculate u(x) for x= -1/2 , -1/4,0,1/4,1/2
        // INDSÆT RUTINER TIL AT BRUGE RESULTATET AF xVec:

        // SLUT RUTINER TIL AT BRUGE RESULTATET AF xVec:
        util::print(xVec);
    }
    // JOHN SLUT
    return 0;
}