function solver(
    Nt,
    dt,
    L,
    N,
    phi,
    x,
    x_dot,
    parameter_phi,
    parameter_impurity,
    parameter_V,
    phi_fct,
    impurity_fct,
    V_fct,
    times_save,
)
    V = [zero(phi) for _ = 1:length(size(phi))+1]
    xarray = [i for i in range(-L / 2.0, L / 2.0, length = N)]
    update_V = ((V, x, t) -> V_fct(V, x, t, xarray, parameter_V...))
    phi_step = ((phi, V, t) -> phi_fct(phi, V, t, dt, parameter_phi...))
    x_step = (
        (phi, V, x, xdot, t) ->
            impurity_fct(phi, V, x, xdot, t, dt, parameter_impurity...)
    )
    phi_r = [zero(phi) for _ = 1:length(times_save)]
    x_r = [zero(x) for _ = 1:length(times_save)]
    x_dot_r = [zero(x_dot) for _ = 1:length(times_save)]
    j = 1
    for i = 1:Nt
        update_V(V, x, i * dt)
        x_step(phi, V, x, x_dot, i * dt)
        phi_step(phi, V, i * dt)
        #push!(o,observables(phi, x, xdot, i*dt)) # this is slow think about it!
        #push!(o,observables(phi, x , xdot,t))
        if i in times_save
            phi_r[j] = deepcopy(phi)
            x_r[j] = deepcopy(x)
            x_dot_r[j] = deepcopy(x_dot)
            j += 1
        end
    end
    return phi_r, x_r, x_dot_r
end

if abspath(PROGRAM_FILE) == @__FILE__
    include("helper.jl")
    include("SGPE.jl")
    include("NewtonImp.jl")
    using BenchmarkTools
    phi = ones(100, 100) + 0im * ones(100, 100)
    V = [
        ones(100, 100) + 0im * ones(100, 100),
        ones(100, 100) + 0im * ones(100, 100),
        ones(100, 100) + 0im * ones(100, 100),
    ]
    xarray = [i for i in range(-5 / 2.0, 5 / 2.0, length = 100)]
    stencil_bc = tuple_stencil_bc(phi)
    dt = 0.1
    lambda = 10.0
    m_b = 1.0
    mu = 0.1
    g = 1.0
    xstep = xarray[2] - xarray[1]
    T = 0.0
    M = 1.0
    re = 1.0
    parameter_phi = [lambda, m_b, mu, g, stencil_bc, xstep, T]
    parameter_impurity = [M, lambda, xstep]
    parameter_V = [re]
    N = 100
    L = 5.0
    x = [0.0, 0.0]
    x_dot = [0.0, 0.0]
end
