How to add, subtract, and take linear combinations of ODE solutions

Hi, I’m solving differential equations with different sets of initial conditions, and I’m trying to figure out how to take different linear combinations of the ODE solutions (added, subtracted, multiplied by different scalar factors, etc.).

For example, with solutions obtained with code like:

prob1 = ODEProblem(DiffEquation!, InitConditionsArray1, tspan, ODEParams1)
prob2 = ODEProblem(DiffEquation!, InitConditionsArray2, tspan, ODEParams2)
soln1 = DifferentialEquations.solve(prob1, Rodas5(), dense=true, reltol=1e-12, abstol=1e-12, maxiters = 1e7, progress=true)
soln2 = DifferentialEquations.solve(prob2, Rodas5(), dense=true, reltol=1e-12, abstol=1e-12, maxiters = 1e7, progress=true)

suppose I want to make a combination solution: solnCombo = 1.5*soln1 - 0.7*soln2.
Is there an easy way to produce this combination as a new variable of the same ODESolution type?

Some notes:
(1) The initial conditions are arrays (over spatial variable x), so these solutions are actually 2-dimensional solutions, functions over a rectangular range of values in (t,x).

(2) It’s probably easy to save the linear combination in an array by just adding/subtracting each data point in the ODESolution's. But, those solutions also have interpolation information in them, so I’m wondering if it’s possible to take linear combinations of the whole solutions.

Thanks for any info!

You can extract arrays of state trajectories using either soln1(t) or Array(soln1) and those can be added if the time steps are the same. You can enforce this using (I think) the keyword saveat. I assume you are talking about the state trajectories when you’re saying you want to form linear combinations of solutions?

Thanks for the comments baggepinnen.

Actually I was hoping to keep all the properties (such as the interpolations) of the ODESolution in the linear combinations, not to turn them into simple arrays. For example, my question is about looking for a way to make the combination solution, solnCombo = 1.5*soln1 - 0.7*soln2, be the same thing as I would’ve gotten from this code:

InitConditionsArrayCombo = 1.5*InitConditionsArray1 - 0.7*InitConditionsArray2
probCombo = ODEProblem(DiffEquation!, InitConditionsArrayCombo, tspan, ODEParams)
solnCombo = DifferentialEquations.solve(probCombo, Rodas5(), dense=true)

The issue is that solving the differential equation can take hours each time, but I may need to take hundreds of linear combinations as part of an optimization process. That’s why I’m looking for a way to compute the ODESolutions soln1 and soln2 just a single time, and take linear combinations of them into a new ODESolution after the fact. (Complete with interpolation properties from the solutions.) Is there any way to do that? Thanks.