@Markiyan Hirnyk Yes, making assumptions may help, but this version

value(I5) assuming x0>0;

produces the result I mentioned.

Your version most likely just evaluates the signum factor:
res:=evalc(value(I5)) assuming x0>0;
indets(res,name);
         {x0, infinity}
evalf(eval(res,x0=1));
         Float(infinity)+Float(infinity)*I
  

@torabi I see that you completely ignored my advice given above.

@torabi You have the following lines, but you don't look at the output. You should!

indets(sys,specfunc(diff)); #Semicolon
## and later
indets(newsys,specfunc(diff));

Notice that the order of the highest derivative of f6 in sys is 2, whereas it is 3 for newsys.
You should avoid differentiating any equation that will increase any of the highest orders of f1,f2, ..., f6.

Right click on the link you have to Steen Markvorsen's file. Choose save as. Open it from within Maple. That worked fine for me.

Let me add that you can solve numerically:

ode := diff(Y(x), x$4) = a^4*Y(x);
ics := Y(0) = 0, D(Y)(0) = 0, Y(l) = 0, (D@@2)(Y)(l) = 0;
ode4:=subs(a^4=a4,ode);
resnum:=dsolve(eval({ode4,ics,(D@@2)(Y)(0)=1},l=2),numeric);
resnum(0);
subs(resnum(0),a4);
%^(1/4); #Same as a2 in the symbolic version
plots:-odeplot(resnum,[x,Y(x)],0..2);
#You can get get the next one by giving an approximate solution:
resnum2:=dsolve(eval({ode4,ics,(D@@2)(Y)(0)=1},l=2),numeric,approxsoln=[a4=4^4,Y(x)=sin(Pi*x)]);
resnum2(0);
subs(resnum2(0),a4);
%^(1/4); #Same as a22
plots:-odeplot(resnum2,[x,Y(x)],0..2);

@patient Since your system is somewhat difficult to examine analytically it may be worthwhile to consider the following example, which can be analyzed analytically and which exhibits similar nonuniqueness features.

restart;
ode1:=diff(v(x),x)=-2*u(x)*sqrt(v(x));
ode2:=v(x)*diff(u(x),x)=0;
dsolve({ode1,ode2,v(0)=1,u(0)=1});
simplify([%]) assuming x<1;
odetest({v(x)=0,u(x)=u(x)},[ode1,ode2]);
#Now picking as an example u(x) as a polynomial or order three and making sure it meets u(x)=1 and u(x)=-1 differentiably at x=1 and 2, respectively:
f:=x->a*x^3+b*x^2+c*x+d;
solve({f(1)=1,D(f)(1)=0,f(2)=-1,D(f)(2)=0},{a,b,c,d});
U12:=eval(f(x),%);
##Piecing a solution together:
U:=piecewise(x<1,1,x<2,U12,-1);
V:=piecewise(x<1,(x-1)^2,x<2,0,(x-2)^2);
plot([V,U],x=0..3,thickness=3,legend=[v(x),u(x)]);

odetest({u(x)=U,v(x)=V},[ode1,ode2]);
#########
## The numerical solution must stop just before 1:
res:=dsolve({ode1,ode2,v(0)=1,u(0)=1},numeric);
plots:-odeplot(res,[[x,v(x)],[x,u(x)]],0..3);

@patient 

1. Using events=[[v(x)=0,halt]] will not have any effect because dsolve/numeric sees the singularity before v(x) = 0 is discovered by 'events'.


It doesn't really make sense to continue the solution after the singularity at x0 as v(x)=0 and u(x)=0 because u(x0) > 0 (and not at all negligible). Thus u will not even be continuous at x0.
What does make sense is to continue the solution for v as v(x) = 0 and then u(x) as anything you like as long as it is differentiable at x0. As we can see from SYS, if v(x)=0 on an interval then it doesn't matter what u(x) is.
Then you will have a solution (v(x), u(x)) that is differentiable. v may not be twice differentiable at x0 though.

2. If you accept that v may not be twice differentiable at x0 then clearly you have nonuniqueness of solutions to the initial value problem on the interval, say -2..2.

Finally, the argument so far relies on the numerical evidence that v(x) and diff(v(x),x) both are zero at x0. That should be examined further.

@patient If you continue v(x) after the singularity (where v(x)=0 and diff(v(x),x) = 0) as v(x) = 0, then u(x) is allowed to be anything.
At x=-.32359993 we have
p(-.32359993);
 u(-.32359993) = 2.73502343609418, v(-.32359993) = 1.32767909590572*10^(-17), D(v)(-.32359993) = -3.07963458440924*10^(-8)

and if v(x)=0 after that then

eval(SYS,v(x)=0);
     {0=0}

thus u(x) can be anything! Just continue it as you please, preferably as a differentiable function.

I don't see how you can avoid the singularity that is actually there. Do you just want to continue from the singularity at x=-.32434965 with u(x) = 0, v(x)=0 ?

@farzane What I start by doing is to show that the equation eq is satisfied by both y=y1=yc*f/(2*E) and by y=y2=yc*x^2. Thus if eq has a unique solution for each given x and ya, then there is no solution (ya(x),yc(x)) to the system of integral equations {eq1, eq2}.
Finally, I give graphical evidence that the equation eq has a unique solution for ya and x given.

@Al86 I'm not sufficiently motivated to use more effort on this topic at least for now. Partly because I shall move to my winter quarters in a few days just like the birds.

@MDD Here is a simple example:

restart;
for k from 1 to 10 do is(b<>0) assuming b^k<>0 end do;

You notice that you receive 'true' for k=1,2,4  but FAIL otherwise.

You might ask, how can 'is' be in doubt whether b=0 or not, if e.g.  b^3<>0  ?
Allowing for the fact that in the short session above b could stand for anything, I cannot conceive of a situation where the third power of whatever could be different from zero if the 'whatever' itself is zero (in the proper meaning of zero in the 'whatever' context).

restart;
assume(b^3<>0);
is(b<>0);
    FAIL

## I think that inequalities (i.e. < ,>, not <> ) are handled much better by the system:
restart;
assume(b^3>0);
is(b<>0); #Answer: true
is(b>0); #Answer: FAIL, which is fine in view of e.g. this example:

b:=exp(2*Pi/3*I);
expand(b^3); # =1 and thus > 0.

@Carl Love The color of the the curly bracket is grey in %piecewise and blue in the other two.

@Carl Love Agree. Indeed the printing could be done by a direct use of `print/piecewise` :

`print/piecewise`(``, x+y-2*z = 1, ``, 2*x+y-3*z = 5, ``, -x+y+z = 1);

or by using unevaluation quotes around piecewise:

'piecewise'(``, x+y-2*z = 1, ``, 2*x+y-3*z = 5, ``, -x+y+z = 1);

@Kitonum The simpler

add(add(a[i,j]*f[i,j], j=1..2), i=1..2);

works too.