Have you looked in the help pages associated with dsolve/numeric starting with

?dsolve,numeric

I'm not sure that the following is what you are thinking about, but you could try it:

plots:-animate(plot,[[x^2+y,x*y,x=-1..1]],y=-1..1);
plots:-animate(plot,[[x^2+y,x*y,x=-1..1]],y=-1..1,trace=24);

sol := solve({x+y=a, x-y=b}, {x,y});

subs(sol,x/y);
#Alternatively
eval(x/y,sol);


Adding an element to a list could be done as shown:

c:= [1,2,4,4,-3];
[op(c),Pi];
[op(1..2,c),Pi,op(3..-1,c)];

A new list is created in each case.

Tables are useful in some cases where the size is not determined beforehand:

c:= [1,2,4,4,-3];
T:=table(c);
T[6]:=Pi;
eval(T);
#Inserting an element meant to be between the elements with indices 3 and 4:
T[3.5]:=hh;
sort(op(eval(T)),(x,y)->lhs(x)<lhs(y));
rhs~(%);

Here are two ways:

X:=[seq(i, i=1..8)];

t:=[2,5,6];

#Method 1:

Y:=X:
for i from 1 to nops(t) do Y:=subs(t[i]=NULL,Y) end do;

#Method 2

eval(subs(`=`~(t,nULL),X),nULL=NULL);

If in the line A:=B you intend to make a new matrix called A and having the entries which B has at that time, then you should use A:=Copy(B) instead.

Avoid assigning to theta[eta] at the end, this will implicitly create a table called theta. Call the result something else. The problem is that you are also using theta as the name of the unknown function.

There will be other problems but that one is critical.

Surely you will have to solve numerically, and I suggest replacing the number 26 by a more modest value to begin with (e.g. 6). Also replace infinity in the boundary conditions with a finite but large value like 1000 (or even less).

Something like this will then work

res:=dsolve({ics, ODE1},numeric,abserr=1e-1,maxmesh=8192);

plots:-odeplot(res,0..1000);

You can use coeffs.

The following lengthy version does not use that the unknowns are indexed except when defining X, at which time you decide their ordering.

restart;

z:= 7*x[1]+11*x[3]-10*x[4]+26*x[6];

X:=[seq(x[i],i=1..6)];
R:=coeffs(z,X,'t');
t;
T:=table(`=`~([t],[R]));

#As an alternative to the elementwise operation ~ you can use zip:
#T:=table(zip(`=`,[t],[R]));
for u in {op(X)} minus {t} do T[u]:=0 end do;
eval(T);
eval(X,op(eval(T)));

Here is a version which makes u[i], DI1, and DI2 functions, and thereby avoids subs.

Also I have made use of `assuming` and removed evalf.

restart;
lambda:=-1; alpha:=1/2;  N:=2;  m:=ceil(alpha);

                             
u[0]:=t->0;      #initial conditions
for i from 0 to N do
 DI1:=unapply(simplify(1/GAMMA(m-alpha)*int((x-t)^(m-alpha-1)*diff(u[i](t),[t$m]),t=0..x)),x) assuming x>0;
 DI2:=unapply(DI1(t)-t^(1/2)/Pi^(1/2)*(u[i](t)+1)^2,t);
 u[i+1]:=unapply(u[i](x)+simplify(int(1/GAMMA(alpha)*lambda*(x-t)^(alpha-1)*DI2(t),t=0..x)),x) assuming x>0;
end do;

Matrix([seq([ n , n^2, n!, n^2 -n!], n= 1..5)]);
latex(%);

Quoting from the help page for solve,

"If the output of the solve command is a piecewise-defined expression, then the assuming command can be used to isolate the desired solution(s). If the output is not piecewise-defined, in particular, if the output is constant, assumptions on the independent variables may be ignored. If there are parameters in the input equations, the solve command will use those assumptions in its computations."

Try this instead:

solve( {x^2=4,x>0},x);

By separating variables in the two ode's (ode1 and ode2) you obtain from E = 1, you find the solutions y(x) = sin(arcsinh(x)) and y(x) = - sin(arcsinh(x)).

However, there are infinitely many solutions since y = 1 and y = -1 are constant solutions which are reached in finite "time" x.

The solution y(x) = sin(arcsinh(x)) reaches 1 for x = x1 = sinh(Pi/2), and infinitely many solutions can be constructed by beginning with sin(arcsinh(x)) until x = x1, then let y(x) = 1 for x1 < x < x2, then switch to a decreasing solution obtained from ode2. And when that reaches -1 you can let the solution be -1 on an interval of arbitrary length then switch (if you like) to an increasing solution obtained from ode1, etc.

Solving numerically your solution of ode1 stays at 1 after it reaches it.

You may try the following though.

restart;
E:=(1+x^2)*(diff(y(x),x))^2+(y(x))^2=1;
E1,E2:=solve(E,{diff(y(x),x)});
E1:=simplify(E1) assuming real;
E2:=simplify(E2) assuming real;
res1:=dsolve(E1 union {y(0)=0});
#This solution is not a solution of E2, but is a solution of E:
simplify(eval(E,res1));
plot(rhs(res1),x=0..12);

#Numerical solution:

res1num:=dsolve(E1 union {y(0)=0},numeric);

plots:-odeplot(res1num,[x,y(x)],0..12);

But as I said there are infinitely many solutions to your problem.

@amrramadaneg You can assign a value to Digits (which initially has the value 10):

Digits:=50;

Now the result of the command

evalf(sin(2)+exp(1));

is 3.6275792552847269307563073372644073404595020651479

whereas

evalhf(sin(2)+exp(1));

is unaffected by the assignment to Digits and returns

3.62757925528472658

round(evalhf(Digits)) gives you the number of digits your computer is using when computing with hardware floats.

For my computer here the value is 15.

See ?evalhf

restart;
LK1 := [1,2,3,4,5]:
#Not knowing the values of the other constants (and since you want a numerical solution) I just set them at some simple values:

K2:=1: nbo4:=2: PH2O:=3: Ti:=4:
seq(fsolve(2*x^3+(2*K1*nbo4+3*K2*nbo4*PH2O)*x^2+3*K1*x*nbo4^2*K2*PH2O-3*K1*nbo4^2*Ti*K2*PH2O = 0, x, 0 .. 10^10),K1=LK1);

(the results in this case are: 1.772001873, 2.117366598, 2.302037223, 2.420167073, 2.503066262).

Notice that for this third order polynomial equation solve can be used to find the 3 exact solutions:

restart;

solve(2*x^3+(2*K1*nbo4+3*K2*nbo4*PH2O)*x^2+3*K1*x*nbo4^2*K2*PH2O-3*K1*nbo4^2*Ti*K2*PH2O = 0, x);