NDSolve for non-analytical forms

Question

I came across a problem of coupled differential equations of a non-analytical results. It follows that NDsolve requires an evaluated form of equations to be fed into it for a differential equation to process.But, my functions can't be evaluated analytically. I want to know if there is any possible way to tackle the problem. I have created a sample problem here:

First I create a Complex system just to confirm I will get non-analytical solutions:

Mat[n_, x_, y_] := SparseArray[{Band[{1, 1}, {n, n}] -> {x^3, x I + 5 y^2 + 4, Sqrt[x]},
                   Band[{1, 2}, {n, n}] -> {y^3, Sqrt[x + I x^3 - y^2 + 4], 
                   Sqrt[x - y^2]}, Band[{3, 1}, {n, n}] -> {I x^3, x + y x^2 + 4, Sqrt[x + y^2]}}]

eval[n_, x_, y_] := Eigensystem[Mat[n, x, y]][[1]]

evec[n_, x_, y_] := Eigensystem[Mat[n, x, y]][[2]]

Secondly, I will form some time-dependent functions so as to cover the complexity of my problem:

val = D[Mat[8, x, y], x];
x2[x1_, t_] := x1 + t^2
y2[y1_, t_] := y1 + t
T1[x1_, y1_] := {I Conjugate[#], #} &@(Conjugate[evec[8, x1, y1][[6]]].SparseArray[
                ArrayRules[val] /. {x -> x1, y -> y1}, Dimensions[val]].evec[8, x1, y1][[7]] // N)
T2[t_] := {t, t^2}
T3[x1_, y1_, t_] := T1[x2[x1, t], y2[y1, t]]
T4[x1_, y1_, t_] := Re[eval[8, x2[x1, t], y2[x1, t]][[7]]] - Re[eval[8, x2[x1, t], y2[x1, t]][[8]]]//N

Last, I will list them in the required form of equation and initial conditons for processing:

t0=-5;
eqns[x1_, y1_, t_] := {A1'[t] == (T2[t].T3[x1, y1, t]) A2[t], A2'[t] == A1[t] (T2[t].Conjugate[T3[x1, y1, t]]),
                      i'[t] == T4[x1, y1, t],A1[t0] == 0, A2[t0] == 1, i[t0] == 0}
sol1 = ParametricNDSolve[eqns[x1, y1, t], {A1[t], A2[t]}, {t, t0, 5}, {x1, y1}]

As you see, say, eqns[x1, y1, t] can't be evaluated unless you provide numerical values of all the parameters. How do we solve the equations in that case. I would be grateful for your help.

(Note: this is a sample just to reflect my problem, feel free to make reasonable changes)

Answer 1

Your equations are of the form x' = f(x), x = {x1,...,xn}, flow part f(x) can be defined implicitly, i.e. as a black box outside of NDSolve.

This explicit form:

sol = NDSolve[{x'[t]==-y[t]-x[t]^2,y'[t]==2x[t]-y[t]^3,x[0]==y[0]==1},{x,y},{t,20}]
ParametricPlot[Evaluate[{x[t],y[t]}/.sol],{t,0,20}]

can be converted to:

ClearAll[flow] ;
flow[x_,y_] := flow[x,y] = {-y-x^2,2 x-y^3} ;

ClearAll[fx,fy] ;
fx[arg__?NumericQ] := flow[arg][[1]] ;
fy[arg__?NumericQ] := flow[arg][[2]] ;

sol = NDSolve[{x'[t] == fx[x[t],y[t]] ,y'[t]==fy[x[t],y[t]],x[0]==y[0]==1},{x,y},{t,20}] ;
ParametricPlot[Evaluate[{x[t],y[t]}/.sol],{t,0,20}]

Edit

In your case:

ClearAll[flow] ;
flow[x1_, y1_, t_,A1_,A2_,i_] :=  flow[x1,y1,t,A1,A2,i] = {
    (T2[t].T3[x1, y1, t]) A2,
    A1 (T2[t].Conjugate[T3[x1, y1, t]]),
    T4[x1, y1, t]
} ;

ClearAll[f1,f2,f3] ;
f1[arg__?NumericQ] := flow[arg][[1]] ;
f2[arg__?NumericQ] := flow[arg][[2]] ;
f3[arg__?NumericQ] := flow[arg][[3]] ;

t0=-5;
(* return only i *)
sol = ParametricNDSolveValue[
    {
        A1'[t] == f1[x1,y1,t,A1[t],A2[t],i[t]],
        A2'[t] == f2[x1,y1,t,A1[t],A2[t],i[t]],
        i'[t] == f3[x1,y1,t,A1[t],A2[t],i[t]],
        A1[t0] == 0, A2[t0] == 1, i[t0] == 0
    },
    i,
    {t, t0, 5},
    {x1, y1}
] 
sol[1,1]