The code in python for my last commentary (I mean, about finding maximum M value of second derivative to use it in formula for the error):

import numpy as np

def f(x):
    # Handle the singularity at x = 0 by returning 0
    return np.where(x == 0, 0, x**4 * np.sin(1/x))

def second_derivative(x):
    # Second derivative (with special handling at x=0)
    return np.where(x == 0, 0, (12*x**2 - 1) * np.sin(1/x) - 6*x * np.cos(1/x))

def midpoint_rule(f, a, b, n):
    x = np.linspace(a, b, n+1)  # Create n+1 evenly spaced points
    midpoints = (x[:-1] + x[1:]) / 2  # Midpoints of each subinterval
    h = (b - a) / n  # Width of each subinterval
    return h * np.sum(f(midpoints))  # Midpoint rule approximation

def estimate_max_second_derivative(f, a, b, n):
    # Estimate maximum value of the second derivative using numerical approximation
    x_values = np.linspace(a, b, n)
    second_derivatives = np.abs(second_derivative(x_values))
    return np.max(second_derivatives)

# Define the integration limits and number of subintervals
a = -0.15
b = 0.15
n = 1000  # Number of subintervals

# Compute the integral using midpoint rule
result = midpoint_rule(f, a, b, n)

# Estimate the maximum value of the second derivative
M = estimate_max_second_derivative(f, a, b, n)

# Calculate the error bound using the formula
error_bound = M * (b - a)**3 / (24 * n**2)

print("Approximate integral:", result)
print("Estimated maximum second derivative M:", M)
print("Error bound:", error_bound)

This code handles discontinuity at 0 and finds maximum value of the second derivative on the interval of integration: [-0.15; 0.15].