Restructure scheme generation

[finsberg]

We would like to make it easier to call into other methods. For example within a scheme we might want to make it possible to call the rhs function.

With the current changes we could implement forward euler as follows

import sympy

def forward_euler(
    ode,
    dt,
    printer, 
    rhs,
    name: str = "values",
    remove_unused: bool = False,
) -> list[str]:
    
    eqs = []

    ret = sympy.IndexedBase(name, shape=(rhs.num_return_values,))
    eqs.append(printer(ret, rhs.states +dt * rhs(states=rhs.states, t=rhs.t, parameters=rhs.parameters), use_variable_prefix=True))
    
    return eqs

or we could implement RK - type schemes, e.g

def heun(
    ode,
    dt,
    printer, 
    rhs,
    name: str = "values",
    remove_unused: bool = False,
) -> list[str]:
    
    eqs = []

    left = sympy.IndexedBase("left", shape=(rhs.num_return_values,))
    right = sympy.IndexedBase("left", shape=(rhs.num_return_values,))
    ret = sympy.IndexedBase(name, shape=(rhs.num_return_values,))
    
    eqs.append(printer(left, rhs(states=rhs.states, t=rhs.t, parameters=rhs.parameters), use_variable_prefix=True))
    eqs.append(printer(right, rhs(states=rhs.states + left * dt, t=rhs.t + dt, parameters=rhs.parameters), use_variable_prefix=True))
    eqs.append(printer(ret, rhs.states + 0.5 * dt * (left + right), use_variable_prefix=True))
    
    return eqs

The challenge now is that we need a reliable way to print objects on the left hand side that are of type IndexedBase. In python this is no problem since we can rely on numpy arrays e.g, the forward euler method here should probably translate into something like

def forward_euler(states, t, dt, parameters):
    return states + dt * rhs(states, t, parameters)

However, in C we might need to use some pointer magic, e.g the forward euler method here should probably translate into something like

void forward_euler(const double *__restrict states, const double t, const double dt, const double *__restrict parameters, double *values) {
    // Assuming size is 10
    double values_0[10];
    rhs(states, t, dt, parameters, values_0);
    for (int i = 0; i < 10; i++){
        values[i] = states[i] + dt * values_0[i] *;
    }
}

or something, like that, so this should be handled separately by the different code generators