Give python code for a 'random walk field solver' to compute capacitance.

Here is just an outline code. It does not work as expected.

import numpy as np
import random
import matplotlib.pyplot as plt

class RandomWalkCapacitanceSolver:
    def __init__(self, grid_size, num_walks, plate_distance, plate_area, boundary_conditions):
        self.grid_size = grid_size
        self.num_walks = num_walks
        self.plate_distance = plate_distance
        self.plate_area = plate_area
        self.boundary_conditions = boundary_conditions
        self.potential_field = np.zeros(grid_size)

        # Initialize boundary conditions
        self.potential_field[0, :] = boundary_conditions['top']
        self.potential_field[-1, :] = boundary_conditions['bottom']

        # Initialize the interior points to the average of the boundary conditions
        interior_potential = (boundary_conditions['top'] + boundary_conditions['bottom']) / 2
        self.potential_field[1:-1, :] = interior_potential

    def solve(self):
        """
        Performs random walks to estimate the potential field and calculates the capacitance.
        Returns both the capacitance and the electric field strength.
        """
        walk_counts = np.zeros(self.grid_size)

        for _ in range(self.num_walks):
            self.perform_walk(walk_counts)

        # Calculate the average potential in the interior
        mask = (walk_counts > 0)
        self.potential_field[mask] /= walk_counts[mask]

        # Calculate electric field strength (V/m)
        voltage_difference = self.boundary_conditions['top'] - self.boundary_conditions['bottom']
        electric_field = voltage_difference / self.plate_distance

        # Calculate capacitance using both methods
        theoretical_capacitance = self.calculate_theoretical_capacitance()
        numerical_capacitance = self.calculate_numerical_capacitance()

        return theoretical_capacitance, numerical_capacitance, electric_field

    def perform_walk(self, walk_counts):
        """
        Performs a single random walk starting from a random point in the grid.
        """
        row, col = self.get_random_start_point()
        potential = 0

        while not self.is_at_boundary(row, col):
            walk_counts[row, col] += 1

            # Update potential based on neighbors
            potential = self.update_potential(row, col)
            self.potential_field[row, col] += potential

            # Get next step
            row, col = self.get_next_step(row, col)

    def update_potential(self, row, col):
        """
        Updates the potential at a point based on its neighbors (Laplace equation).
        """
        neighbors = []
        
        # Only add valid neighbors within bounds
        if row > 0:     neighbors.append(self.potential_field[row-1, col])  # Up
        if row < self.grid_size[0] - 1: neighbors.append(self.potential_field[row+1, col])  # Down
        if col > 0:     neighbors.append(self.potential_field[row, col-1])  # Left
        if col < self.grid_size[1] - 1: neighbors.append(self.potential_field[row, col+1])  # Right
        
        # Return the average of the valid neighbors
        return sum(neighbors) / len(neighbors) if neighbors else 0

    def calculate_theoretical_capacitance(self):
        """
        Calculates the theoretical capacitance using the parallel plate formula.
        """
        epsilon_0 = 8.85e-12  # Permittivity of free space in F/m
        return epsilon_0 * self.plate_area / self.plate_distance

    def calculate_numerical_capacitance(self):
        """
        Calculates the numerical capacitance based on the potential field.
        """
        # Calculate the electric field at the plate surface
        surface_field = np.abs(self.potential_field[1] - self.potential_field[0]) / (self.plate_distance / self.grid_size[0])

        # Calculate the charge density (σ = ε₀E)
        epsilon_0 = 8.85e-12
        charge_density = epsilon_0 * surface_field

        # Calculate total charge (Q = σA)
        total_charge = charge_density * self.plate_area

        # Calculate capacitance (C = Q/V)
        voltage_difference = self.boundary_conditions['top'] - self.boundary_conditions['bottom']
        return total_charge / voltage_difference if voltage_difference != 0 else 0

    def get_random_start_point(self):
        """Returns a random starting point in the interior of the grid."""
        row = random.randint(1, self.grid_size[0] - 2)
        col = random.randint(0, self.grid_size[1] - 1)
        return row, col

    def get_next_step(self, row, col):
        """Returns the next valid step in the random walk."""
        while True:
            direction = random.randint(0, 3)
            next_row, next_col = row, col

            if direction == 0:   next_row -= 1  # Up
            elif direction == 1: next_row += 1  # Down
            elif direction == 2: next_col -= 1  # Left
            elif direction == 3: next_col += 1  # Right

            # Handle periodic boundary conditions for left/right
            if next_col < 0:
                next_col = 0
            elif next_col >= self.grid_size[1]:
                next_col = self.grid_size[1] - 1

            # Handle top and bottom boundaries
            if next_row < 0:
                next_row = 0
            elif next_row >= self.grid_size[0]:
                next_row = self.grid_size[0] - 1

            return next_row, next_col

    def is_at_boundary(self, row, col):
        """Checks if the current position is at a boundary (plate)."""
        return row == 0 or row == self.grid_size[0] - 1

    def plot_potential_field(self):
        """Plots the potential field using a heatmap."""
        plt.figure(figsize=(10, 8))
        im = plt.imshow(self.potential_field, cmap='viridis', interpolation='nearest')
        plt.colorbar(im, label='Potential (V)')
        plt.title('Potential Field between Parallel Plates')
        plt.xlabel('Columns')
        plt.ylabel('Rows')
        plt.show()

# Example usage
grid_size = (50, 50)
num_walks = 100000  # Increased for better convergence
plate_distance = 0.01  # 1 cm
plate_area = 1.0    # 1 m²
boundary_conditions = {
    'top': 1.0,     # Top plate: 1V
    'bottom': 0.0   # Bottom plate: 0V
}

solver = RandomWalkCapacitanceSolver(grid_size, num_walks, plate_distance, plate_area, boundary_conditions)
theoretical_cap, numerical_cap, e_field = solver.solve()

print(f"Theoretical capacitance: {theoretical_cap:.2e} F")
print(f"Numerical capacitance: {numerical_cap:.2e} F")
print(f"Electric field strength: {e_field:.2e} V/m")

solver.plot_potential_field()