Równaie poissona kłopot z deklaracją.

Mógłby ktoś mi pomóc i uswiadomić w którym miejcu i w jaki sposób zadeklarowac funkcję v9x,y0=6xy(1-y)-2x^3 i
i warunki brzegowe dla x=0 : alfa 1=1, beta1=0, gamma1=0, a dla x=1 alfa2=1, beta2=0, i gamma2=y(1-y).

Kod rozwiązywania rówananie Piossona 2d mam następujący

#include <cmath>
#include <cstdlib>
#include <ctime>
#include <iostream>
#include <fstream>
#include <vector>
using namespace std;

int L;                        
vector< vector<double> > psi; // potencjał
vector< vector<double> > rho; //gęstośc ładunku- podana
double h;                     // krok h
vector< vector<double> > psiNew;   // potencjał po każdym nowym kroku
int steps;                    // liczba iteracji
double accuracy;              // dokładność


void initialize() {

    // tworzenie macierzy
    int N = L + 2;            // punkty wewnętrzne i dwa warunki brzegowe
    psi = psiNew = rho
      = vector< vector<double> >(N, vector<double>(N, 0.0));

    h = 1 / double(L + 1);    // założenie wymiaru x,y=1
    double q = 10;            // point charge
    int i = L / 2;            // środek siatki
    rho[i][i] = q / (h * h);  // gęstośc ładunku

    steps = 0;
}

void Jacobi() {

    // Metoda Jacobiego dla pojedynczej iteracji
    for (int i = 1; i <= L; i++)
    for (int j = 1; j <= L; j++)
        psiNew[i][j] = 0.25 * (psi[i - 1][j] + psi[i + 1][j] +
                               psi[i][j - 1] + psi[i][j + 1] +
                               h * h * rho[i][j]);
}

double relativeError() {

    double error = 0;         // średni błąd 
    int n = 0;                // liczba niezerowych różnic

    for (int i = 1; i <= L; i++)
    for (int j = 1; j <= L; j++) {
        if (psiNew[i][j] != 0)
        if (psiNew[i][j] != psi[i][j]) {
            error += abs(1 - psi[i][j] / psiNew[i][j]);
            ++n;
        }
    }
    if (n != 0)
        error /= n;

    return error;
}

void GaussSeidel() {

    // copy psi to psiNew
    for (int i = 1; i <= L; i++)
    for (int j = 1; j <= L; j++)
        psiNew[i][j] = psi[i][j];

    // Gauss-Seidel update
    for (int i = 1; i <= L; i++)
    for (int j = 1; j <= L; j++)
        psiNew[i][j] = 0.25 * (psiNew[i - 1][j] + psiNew[i + 1][j] +
                               psiNew[i][j - 1] + psiNew[i][j + 1] +
                               h * h * rho[i][j]);
}



void iterate(void (*method)()) {

    while (true) {
        method();
        ++steps;
        double error = relativeError();
        if (error < accuracy)
            break;
        vector< vector<double> > swap = psi;
        psi = psiNew;
        psiNew = swap;
    }
    cout << " :Liczba kroków = " << steps << endl;
}

int main() {

    cout << " Iteracyjne rozwiązanie równania Poissona\n"
         << " ----------------------------------------\n";
    cout << " Wprowadź liczbę punktów x i y  : ";
    cin >> L;

    initialize();

    cout << " Wprowadź dokładniość ";
    cin >> accuracy;

    clock_t t0 = clock();
    cout << " Wybierz metodę Jacobiego 1 lub Gaussa siedla 2 ";
    int choice;
    cin >> choice;
    switch (choice) {
    case 1:
        iterate(Jacobi);
        break;
    case 2:
        iterate(GaussSeidel);
        break;
   
    default:
        cout << " Jacobi: " << endl;
        iterate(Jacobi);
        cout << " Gauss-Seidel: " << endl;
        initialize();
        iterate(GaussSeidel);
       
        initialize();
       
        break;
    }

    clock_t t1 = clock();
    cout << " CPU time = " << double(t1 - t0) / CLOCKS_PER_SEC
         << " sec" << endl;

    
    cout << " Potenciał w folderze  poisson.data" << endl;
    ofstream file("poisson.data");
    for (int i = 0; i < L + 2; i++) {
        double x = i * h;
        for (int j = 0; j < L + 2; j++) {
            double y = j * h;
            file << x << '\t' << y << '\t' << psi[i][j] << '\n';
        }
        file << '\n';
    }
    file.close();
}