def METHOD_NAME(mesh, layers, bias, init_method="normal"):
    init_method = init_method.lower()
    assert init_method in ["normal", "uniform"]
    weights = []
    for i in range(len(layers) - 1):
        R_layer, vector_layer = max(layers[i + 1], layers[i]), min(layers[i + 1], layers[i])
        R = VectorFunctionSpace(mesh, "R", 0, dim=R_layer)
        ws = []
        weight = {}
        for _ in range(vector_layer):
            W = Function(R)
            if init_method == "uniform":
                W.vector()[:] = random(R.dim())
            elif init_method == "normal":
                W.vector()[:] = np.sqrt(2 / vector_layer) * randn(R.dim())
            ws.append(W)
        weight["coefficient"] = ws
        if bias[i]:
            R = VectorFunctionSpace(mesh, "R", 0, dim=layers[i + 1])
            weight["bias"] = Function(R)
        weights.append(weight)
    return weights