How to remove this numerical artifact?

How to solve this (perhaps a little more complicated than necessary) with the tools of python scipy.integrate I demonstrated in How to numerically set up to solve this differential equation?

---

If you want to stay with the simplicity of a one-stage method, expand the step as $f(t+s)=f(t)+h(s)$ where $t$ is constant and $s$ the variable, so that $$ εh'(s)=εf'(t+s)=g(t+s)-f(t)^4-4f(t)^3h(s)-6f(t)^2h(s)^2-... $$ The factor linear in $h$ can be moved to and integrated into the left side by an exponential integrating factor. The remaining terms are quadratic or of higher degree in $h(Δt)\simΔt$ and thus do not influence the order of the resulting exponential-Euler method. \begin{align} ε\left(e^{4f(t)^3s/ε}h(s)\right)'&=e^{4f(t)^3s/ε}\left(g(t+s)-f(t)^4-6f(t)^2h(s)^2-...\right) \\ \implies h(Δt)&\approx h(0)e^{-4f(t)^3Δt/ε}+\frac{1-e^{-4f(t)^3Δt/ε}}{4f(t)^3}\left(g(t)-f(t)^4\right) \\ \implies f(t+Δt)&\approx f(t)+\frac{1-e^{-4f(t)^3Δt/ε}}{4f(t)^3}\left(g(t)-f(t)^4\right) \end{align}

This can be implemented as

eps = 0.03

def step(t,f,dt):
# exponential Euler step
    g = max(0,np.sin(t))
    f3 = 4*f**3;
    ef = np.exp(-f3*dt/eps)
    return f + (1-ef)/f3*(g - f**4)

# plot the equilibrium curve f(t)**4 = max(0,sin(t))
x = np.linspace(0,np.pi, 150);
plt.plot(x,np.sin(x)**0.25,c="lightgray",lw=5)
plt.plot(2*np.pi-x,0*x,c="lightgray",lw=5)

for N in [500, 100, 50]:
    a0, a1 = 0, eps/2
    t = np.linspace(0,2*np.pi,N+1)
    dt = t[1]-t[0];
    while abs(a0-a1)>1e-6:
        # Aitken delta-squared method to accelerate the fixed-point iteration
        f = a0 = a1;
        for n in range(N): f = step(t[n],f,dt);
        a1 = f;
        if abs(a1-a0) < 1e-12: break
        for n in range(N): f = step(t[n],f,dt);
        a2 = f;
        a1 = a0 - (a1-a0)**2/(a2+a0-2*a1)
    # produce the function table for the numerical solution
    f = np.zeros_like(t)
    f[0] = a1;
    for n in range(N): f[n+1] = step(t[n],f[n],dt);
    plt.plot(t,f,"-o", lw=2, ms=2+200.0/N, label="N=%4d"%N)

plt.grid(); plt.legend(); plt.show()

and gives the plot

showing stability even for $N=50$. The errors for smaller $N$ look more chaotic due to the higher non-linearity of the method.

If you are not told to do it all by yourself, I would suggest you to use the powerful scipy package (specially the integrate subpackage) which exposes many useful objects and methods to solve ODE.

import numpy as np
from scipy import integrate
import matplotlib.pyplot as plt

First define your model:

def model(t, y, c=0.03):
    return (np.max([np.sin(t), 0]) - y**4)/c

Then choose and instantiate the ODE Solver of your choice (here I have chosen BDF solver):

t0 = 0
tmax = 10
y0 = np.array([0.35]) # You should compute the boundary condition more rigorously
ode = integrate.BDF(model, t0, y0, tmax)

The new API of ODE Solver allows user to control integration step by step:

t, y = [], []
while ode.status == 'running':
    ode.step() # Perform one integration step
    # You can perform all desired check here...
    # Object contains all information about the step performed and the current state! 
    t.append(ode.t)
    y.append(ode.y)
ode.status # finished

Notice the old API is still present, but gives less control on the integration process:

t2 = np.linspace(0, tmax, 100)
sol = integrate.odeint(model, y0, t2, tfirst=True)

And now requires the switch tfirst set to true because scipy swapped variable positions in model signature when creating the new API.

Both result are compliant and seems to converge for the given setup:

fig, axe = plt.subplots()
axe.plot(t, y, label="BDF")
axe.plot(t2, sol, '+', label="odeint")
axe.set_title(r"ODE: $\frac{d f}{d\theta} = \frac{1}{c}(\max(\sin\theta, 0) - f^4)$")
axe.set_xlabel("$t$")
axe.set_ylabel("$y(t)$")
axe.set_ylim([0, 1.2])
axe.legend()
axe.grid()

Solving ODE numerically is about choosing a suitable integration method (stable, convergent) and well setup the parameters.

I have observed that RK45 also performs well for this problem and requires less steps than BDF for your setup. Up to you to choose the Solver which suits you best.