One order method

import numpy as np
from scipy.optimize import root
import plotly.graph_objects as go
from plotly.subplots import make_subplots
import plotly.io as pio
pio.templates.default = "seaborn"

class curtiss_model:

    def __init__(self, k):
        self.k = k

    def fcn(self, t, y) :
        k = self.k
        y_dot = k * (np.cos(t) - y)
        return y_dot

    def sol(self, yini, t0, t):
        k = self.k
        c0 = (yini - (k/(k*k + 1)) * (k*np.cos(t0) + np.sin(t0))) * np.exp(-k*t0)
        y = (k/(k*k + 1)) * (k*np.cos(t) + np.sin(t)) +  c0 * np.exp(-k*t)
        return y

yini = 2.
tini = 0.
tend = 2.
k = 10.

cm = curtiss_model(k)
    
texa = np.linspace(tini, tend, 1000)
uexa = cm.sol(yini, tini, texa)

fig = go.Figure()
fig.add_trace(go.Scatter(x=texa, y=uexa, name='f(x)', marker_maxdisplayed=200))

#create slider
steps = []
for k_i in range(10, 101, 10):
    cm = curtiss_model(k_i)
    step = dict(method="update", label = f"{k_i}", args=[{"y": [cm.sol(yini, tini, texa)]}])
    steps.append(step)
sliders = [dict(currentvalue={'prefix': 'k = '}, steps=steps, pad={'t': 50})]

fig.update_layout(sliders=sliders, title = 'Exact solution of the Curtiss and Hirschfelder equation', 
                  xaxis_title="t", yaxis_title="y",  height=500, modebar=dict(orientation='v'))

fig.show()

class ode_result:
    def __init__(self, y, t):
        self.y = y
        self.t = t

def forward_euler(tini, tend, nt, yini, fcn):

    dt = (tend-tini) / (nt-1)
    t = np.linspace(tini, tend, nt)

    yini = np.atleast_1d(yini)
    neq = yini.size

    y = np.zeros((neq, nt))
    y[:,0] = yini

    for it, tn  in enumerate(t[:-1]):
        yn = y[:,it]
        y[:,it+1] = yn + dt*np.atleast_1d(fcn(tn, yn))

    return ode_result(y, t)

def sol_and_error_curtiss(yini, tini, tend, k, nt, method):
    
    if method == forward_euler:
        str_method = "Forward Euler"
    elif method == backward_euler:
        str_method = "Backward Euler"
    else:
        print("Unknown method")
        return
    
    cm = curtiss_model(k)
    fcn = cm.fcn
    
    texa = np.linspace(tini, tend, 500)
    yexa = cm.sol(yini, tini, texa)

    nt_max= nt[-1]
    dt_max = (tend-tini)/(nt_max-1)
    sol = method(tini, tend, nt_max, yini, fcn)
    err = np.abs(sol.y[0]-cm.sol(yini, tini, sol.t))[1:]
    
    # first plot
    fig = make_subplots(rows=2, cols=1, subplot_titles=(str_method, "Gloabl error"), vertical_spacing=0.1)
    color = fig.layout.template.layout.colorway
    fig.add_trace(go.Scatter(x=texa, y=yexa, name='Exact sol.', line_color='grey', legendgroup = 'sol'), 
                  row=1, col=1)
    fig.add_trace(go.Scatter(x=sol.t, y=sol.y[0], mode='markers+lines', line_dash='dot', 
                             line_color=color[0], name='Numercal sol. ', legendgroup = 'sol'), row=1, col=1)
    fig.add_trace(go.Scatter(x=sol.t[1:], y=err, mode='markers+lines', line_dash='dot',
                             line_color=color[1], name='error', legendgroup = 'err'), row=2, col=1)
    
    # create slider
    steps = []
    for i, nt_i in enumerate(nt):
        sol = method(tini, tend, nt_i, yini, fcn)
        err = np.abs(sol.y[0]-cm.sol(yini, tini, sol.t))[1:]
        dt = (tend-tini)/(nt_i-1)
        step = dict(method="update", label = f"{nt_i}", args=[{"x": [texa, sol.t, sol.t[1:]], 
                                                               "y": [yexa, sol.y[0], err]},
                {"title": f"Solution and error for k={int(k)}, dt={dt:.4e} and k.dt={k*dt:.2f}"}])
        steps.append(step)

    sliders = [dict(active=len(nt)-1, currentvalue={"prefix": "nt : "}, steps=steps)]

    fig.update_yaxes(exponentformat='e', col=1)
    legend = dict(tracegroupgap=300, groupclick="toggleitem", x=0.8, bgcolor='rgba(0,0,0,0)')
    fig.update_layout(height=800, sliders=sliders, legend=legend, modebar=dict(orientation='v'),
                      title=f"Solution and error for k={int(k)}, dt={dt_max:.4e} and k.dt={k*dt:.2f}")

    fig.show()

def backward_euler(tini, tend, nt, yini, fcn):

    dt = (tend-tini) / (nt-1)
    t = np.linspace(tini, tend, nt)

    yini = np.atleast_1d(yini)
    neq = yini.size

    y = np.zeros((neq, nt), order='F')
    y[:,0] = yini

    def g(uip1, *args):
        uip, tip1 = args
        return uip1 - uip - dt*np.array(fcn(tip1, uip1))

    for it, tn  in enumerate(t[:-1]):
        yn = y[:,it]
        y0 = yn + dt*np.atleast_1d(fcn(tn, yn))
        # solve y[:,it+1] - y[:,it] - dt * fcn(tini + (it+1)*dt, y[:,it+1]) = 0
        sol = root(g, y0, (yn, tn+dt))
        y[:,it+1] = sol.x

    return ode_result(y, t)