Here is a simple version.

I have assumed that the system itself doesn't contain parameters. Therefore I have kept it outside the procedure.

As it has been made the procedure returns a sequence of 4 numerical procedures.

sys:=diff(f1(t),t)=piecewise(t<=4,3*exp(-2*t),-1),diff(f2(t),t)=sin(t),diff(g1(t),t)=f1(t),diff(g2(t),t)=f2(t);

p:=proc(a1::realcons,a2::realcons)
local init,F;
init:= f1(0)=a1, f2(0)=a2,g1(0)=0,g2(0)=0;
F:=dsolve({sys,init}, numeric,output=listprocedure);
op(subs(F,[f1(t),f2(t),g1(t),g2(t)]));
end proc:
F1,F2,G1,G2:=p(1,1);
G2(5.5)/G1(5.5);
F1,F2,G1,G2:=p(Pi,1);
G2(5.5)/G1(5.5);

plots:-animate(plot,['p(a,1)[1]' ,0..6],a=.5..3);

Perhaps you would prefer this other version, which returns the listprocedure:

p2:=proc(a1::realcons,a2::realcons)
local init,F;
init:= f1(0)=a1, f2(0)=a2,g1(0)=0,g2(0)=0;
dsolve({sys,init}, numeric,output=listprocedure);
end proc:

p2(Pi,1);
F1,F2,G1,G2:=op(subs(%,[f1(t),f2(t),g1(t),g2(t)]));
 

Preben Alsholm
 

I don't see the need for aold. Does  this accomplish what you want?

restart;
anew:=x-> x^2:
for i from 1 to 5 do
anew:=unapply(anew(x)+x,x);
end do:
anew(x);

Preben Alsholm

This works.

restart;
sys:=diff(f1(t),t)=piecewise(t<=4,3*exp(-2*t),-1),
diff(f2(t),t)=sin(t);
init:= f1(0)=1, f2(0)=1:
F:=dsolve({sys,init}, {f1(t), f2(t)}, numeric,output=listprocedure);
F1:=subs(F,f1(t));
F2:=subs(F,f2(t));
int(F2(t),t=0..5.5)/int(F1(t),t=0..5.5);
 

I guess you insist on doing it numerically, but in this case it is certainly not necessary.

Preben Alsholm

Is this what you want?

with(plots):
a := animate(spacecurve, [[cos(t), sin(t), t], t = 0 .. A], A = 0 .. 2*Pi, color = red, axes = boxed, labels = [x, y, z]):
b := animate(spacecurve, [[piecewise(t<Pi,undefined,sin(t)),cos(t), t], t = 0 .. A], A = 0 .. 2*Pi, color = blue):
display(a, b);
 

Preben Alsholm

Depending on which one of n or p you want to eliminate in favor of sigma, you could use solve in combination with subs.

Here I'm eliminating n:

solve(Variance(X3)=sigma^2, {n}):
j3:=subs(%,k3);
 

Preben Alsholm

This method seems to work well.

I have minimized the residuals of the first 10 equations and treated the last two as equality constraints.

This gives an answer, which makes it possible for fsolve to succeed.

res2:=LSSolve(Res[1..10],{eq11,eq12,x1<=1,x2<=1,x3<=1,x4<=1,x5<=1,x6<=1,x7<=1,x8<=1,x9<=1,x10<=1,T>=Ti},
assume=nonnegative,iterationlimit=5000,feasibilitytolerance=.001,optimalitytolerance=0.001);
evalf(eval([eq1, eq2,eq3,eq4,eq5,eq6,eq7,eq8,eq9,eq10,eq11,eq12], res2[2]));
fsolve({eq1, eq2,eq3,eq4,eq5,eq6,eq7,eq8,eq9,eq10,eq11,eq12}, convert(res2[2],set));

{T = 3410.570845, n = 9.711378242, x1 = 0.03121886702, x10 = 0.01450927764,

  x2 = 0.1420646247, x3 = 0.6380365355, x4 = 0.009296381139,

  x5 = 0.07175312883, x6 = 0.04147335723, x7 = 0.01776536923,

  x8 = 0.006835808416, x9 = 0.02704665037}
 

Preben Alsholm

Here is my code after having executed your assignments of equations etc.

(lhs means 'left hand side: see ?lhs)

I have multiplied the last equation by 10^(-2) to reduce its weight in the optimization. This was the least weigt I had success with.

Preben Alsholm

with(Optimization):
Res:=map(lhs-rhs,[eq1, eq2,eq3,eq4,eq5,eq6,eq7,eq8,eq9,eq10,eq11,eq12*10^(-2)]):
#res:=LSSolve(Res,{x1<=1,x2<=1,x3<=1,x4<=1,x5<=1,x6<=1,x7<=1,x8<=1,x9<=1,x10<=1,T>=Ti},assume=nonnegative,iterationlimit=5000);
res:=LSSolve(Res,{x1<=1,x2<=1,x3<=1,x4<=1,x5<=1,x6<=1,x7<=1,x8<=1,x9<=1,x10<=1,T>=Ti},assume=nonnegative,iterationlimit=5000,optimalitytolerance=0.0001);
[0.00517923021408238112, [T = 2782.09143970567812, n = 10.1585637591277820,

  x1 = 0.00872658980230692216, x10 = 0.0748602181855385214,

  x2 = 0.0766666505398920090, x3 = 0.579521666596851892,

  x4 = 0.0000115210794556045040, x5 = 0.0901704459988488805,

  x6 = 0.0852556575745576917, x7 = 0.00577381343641316520,

  x8 = 0.00427475756541988076, x9 = 0.0642723480725257840]]
evalf(eval([eq1, eq2,eq3,eq4,eq5,eq6,eq7,eq8,eq9,eq10,eq11,eq12], res[2]));
[                                                              
[1. = 1.004651844, 4. = 4.001364580, 3.333333334 = 3.329171039,

  12.53333334 = 12.53468790, 0.9895336689 = 1., 0.02234128225 = 0.03162447862,

  0.0001704541989 = 0.02341381146, 0.001029281057 = 0.06427234807,

  0.0002380276970 = 0.07486021819, 0.07729879942 = 0.07666665054,

                                                5                 5]
  0.01208325906 = 0.008726589802, 6.693006871 10  = 6.693006871 10 ]
 

The underscore is often used by Maple when it has to come up with a global name, here a summation index.

In finding the general solution to a differential equation Maple creates names like _C1, _C2, etc.

I think the underscore is meant as a signal that this is a global variable invented by Maple, and don't use that kind of name yourself.

That the names are global implies that you can do e.g.

subs(_t=k,%);
                      infinity                        
                       -----                          
                        \     /   (1 - rho)          k\
                         )    |  k          q (1 - q) |
                        /     |- ---------------------|
                       -----  \        -1 + rho       /
                       k = 0                          
Preben Alsholm

C:=x->a(2*x^2+(V/x))+b(8*x+(V/x^2));
diff(C(x),x);
#The derivative as a function:
D(C);

# But maybe what you really mean is:

C:=x->a*(2*x^2+V/x)+b*(8*x+V/x^2);

Same syntax, though.

Preben Alsholm
 

You may try LSSolve in the Optimization package.

with(Optimization):
Res:=map(lhs-rhs,[eq1, eq2,eq3,eq4,eq5,eq6,eq7,eq8,eq9,eq10,eq11,eq12]):
res:=LSSolve(Res,{x1<=1,x2<=1,x3<=1,x4<=1,x5<=1,x6<=1,x7<=1,x8<=1,x9<=1,x10<=1,T>=Ti},assume=nonnegative,iterationlimit=1000);
Warning, limiting number of major iterations has been reached
  [199.025377293247600, [T = 2913.65979012579737, n = 1.89800495418788406,

    x1 = 0.419279311450406810, x10 = 0.489249736215187026,

    x2 = 0.403106848317799415, x3 = 0.985749025693950176,

    x4 = 0.330160746572949249, x5 = 0.242974425918598997,

    x6 = 0.189383244239871934, x7 = 0.372912480520804446,

    x8 = 0.353892475782923488, x9 = 0.393349253160286405]]

evalf(eval([eq1, eq2,eq3,eq4,eq5,eq6,eq7,eq8,eq9,eq10,eq11,eq12], res[2])):

fsolve({eq1, eq2,eq3,eq4,eq5,eq6,eq7,eq8,eq9,eq10,eq11,eq12}, convert(res[2],set));

The results of the last two lines show that you still don't have a solution to your equations, but at least the right and left hand sides are of the same order of magnitude.

Preben Alsholm

Here is a way to do it.

ode:=diff(y(x),x)=y(x);
sol1:=dsolve({ode,y(0)=1},numeric);
sol1(2);
#In this case op(2,sol1(x)):
Y1 := x-> rhs(op(2, sol1(x)));
Y1(.1234);
evalf(Int(Y1, 0 .. 1));
 

# or instead

sol2:=dsolve({ode,y(0)=1},numeric,output=listprocedure);
Y2:=subs(sol2,y(x));
evalf(Int(Y2, 0 .. 1));
 

My preferred method is the latter.

Preben Alsholm

The following also works, but is a lot slower than the method given by Robert Israel.

evalf(D[2](uu)(2,1));
plot(D[2](uu)(x,2),x=0..5,numpoints=5,adaptive =false);
evalf(Int(D[2](uu)(x,2),x=0..5,method=_Dexp,epsilon=1e-3));

Preben Alsholm
 

Generally speaking, I think this message is really bad news.

I guess you don't have a backup.

It is a good idea to save your document without output.  Remove output before saving.

Hopefully, someone can help you. But my experience tells me that there is not much to do about it.

Preben Alsholm

There is a very good reason why Maple cannot solve for the leading derivatives although they appear linearly in the system RedSys.

Consider the linear system with the unknowns being the leading derivatives, X16_1,X21_1,X63_2,X74_3,X97_2.

Then the coefficient matrix has only got rank 3. Had it been 5, there wouldn't have been a problem.

SYS:=subs(diff(X74(t),t,t,t)=X74_3,diff(X74(t),t,t)=X74_2,diff(X74(t),t)=X74_1,X74(t)=X74,diff(X97(t),t,t)=X97_2,diff(X97(t),t)=X97_1,X97(t)=X97,diff(X63(t),t,t)=X63_2,diff(X63(t),t)=X63_1,X63(t)=X63,diff(X21(t),t)=X21_1,X21(t)=X21,diff(X16(t),t)=X16_1,X16(t)=X16,RedSys):
L:=[op(1..4,SYS),lhs(SYS[5])=rhs5]:
A,b:=LinearAlgebra:-GenerateMatrix(L,[X16_1,X21_1,X63_2,X74_3,X97_2]):
LinearAlgebra:-Rank(A);
LinearAlgebra:-GaussianElimination(A);
LinearAlgebra:-Rank(<A|b>);

Preben Alsholm

Besides maybe other problems there is a problem with the initial conditions:

convert(RedSys,D):
eval(%,t=0):
res:=eval(%,incons1);
   [                                                                 
   [0 = 0.002, -800.0100000 @@(D, 2)(X97)(0) = 0, 0 = 0, 4.106328774 10  

                                                    
     @@(D, 2)(X63)(0) - 6.056424308 10   @@(D, 3)(X74)(0)

                                                      ]
      - 0.000002204071969 @@(D, 2)(X74)(0) = 0, 0. = 0]

Preben Alsholm