@maxwell You will probably find a numerical solution easier to find, but of course that means that all parameters must be concrete and that you won't have a formula.
Starting where I left off above i.e. with constants defined (and with sol3A found for comparison only) try this:

res3:=dsolve({deq3, v3(0) = v__0}, numeric);
plots:-odeplot(res3,[t,v3(t)],0..10);
plots:-odeplot(res3,[t,v3(t)-sol3A],0..10);

Here is the last plot. It represents the difference between the exact solution (sol3A) and the numerical solution:

You can make that error smaller if desired by decreasing the default values of abserr and relerr, which for method=rkf45 are abserr=1e-7 and relerr=1e-6. (See the help page ?dsolve/numeric/rkf45)

Digits:=15:
res3a:=dsolve({deq3, v3(0) = v__0}, numeric,abserr=1e-12,relerr=1e-12);
plots:-odeplot(res3a,[t,v3(t)-sol3A],0..10);

@nm You could try the functions initially known to Maple.
Examples: Replace foo by sin, exp, LambertW, or whatever you can come up with in
showstat(`evalf/foo`);

You can define your own procedure `evalf/myfunc` for your function myfunc.
See ?evalf  (details)

@ I added a second numerical version to my answer (at the bottom). It defines peano as before (adding option remember though). Then it defines two very simple procedures, which just call peano. peano1(...)  returns peano(...)[1] and peano2(...)  returns peano(...)[2].
Notice that if the second argument (depth, here 5) is fixed then peano(t,depth) is only a function of one variable, but peano returns two results. Thus the need for peano1 and peano2. But you cannot expect those two to have common zeros.

Here is a little history:
The bug doesn't appear in Maple 8 nor in Maple 12. (The representations of the solution are different in the two, but they are equivalent).
The bug appears in Maple 15, 16, 17, 18.
The bug doesn't appear in Maple 2015.2 nor in Maple 2016.2. (sin and cos are used for the imaginary roots as they also are in Maple 12).
Whether it is present in Maple 13 or 14 I don't know.
### The following does work in Maple 18 (with a,b,c unassigned):

ode2:=a*diff(y(x), x$4)+b*diff(y(x),x$2)+c*y(x) = 0;
sol2:=dsolve(ode2);
eval(sol2,{a=26,b=-50,c=-2});
simplify(%);

@DuncanNox You wrote:
"What is a large number of n (or dimension of system) for solving mentioned equation in maple? I mean how many steps is possible to do in this metod using todays PC."
I don't know. But you could try with higher values of n.

@DuncanNox We prefer text so we don't have to type.
1. You have a syntax error: It should be s[i]:=subs(Soln,y(x));  since the right hand sides in Soln are procedures.
2. You forgot known=s[i-1].
###
But here is the approach that I called simplest and best:

restart;
n:=50:
sys:=seq(diff(y[i](x),x)=6*x-2*y[i-1](x),i=1..n):
y[0]:=x->3;
ics:=seq(y[i](0)=3,i=1..n):
res:=dsolve({sys,ics},numeric);
plots:-odeplot(res,[seq([x,y[i](x)],i=0..n)],0..3);
#Animation shown below:
plots:-display(seq(plots:-odeplot(res,[x,y[i](x)],0..3,title=i),i=0..n),insequence,labels=[x,y]);
### Symbolic solution:
sol:=dsolve({sys,ics}):
nops(indets(sol,specfunc(Int)));# Unevaluated Int present
sol1:=value(sol): #No more!
nops(indets(sol1,specfunc(Int)));
Y:=subs(sol1,[seq(y[i](x),i=1..n)]):
plot(Y,x=0..3);
plots:-display(seq(plot(Y[i],x=0..3),i=1..n),insequence);

Animating the numerical solution:

@DuncanNox
You wrote: " that would work, while there only two ODE. Bud this is only a part of my problem. A want to take first solution to second ODE, second solution do third ODE, atc. many times. "

dsolve can solve large systems, so what is the problem? Could you give us a concrete example?

@vv You wrote " I don't see why including elif in the definition would be useful, except that g(x) is now undefined for a symbolic x. "
My point was only that IF you want to use is then you must include the elif part. The undefined might as well have been FAIL or 345.

@vv As I see it this is not a bug. Quoting the relevant section from the help page on this in full. (It talks about is(x1,prop1) ).

"The is routine determines whether x1 satisfies the property prop1. It returns true, false, or FAIL. The is routine returns true if all possible values of x1 satisfy the property prop1. The is routine returns false if any possible value of x1 does not satisfy the property prop1. The is function returns FAIL if it cannot determine whether the property is always satisfied. This is a result of insufficient information or an inability to compute the logical derivation."

Note: I have corrected the initial version of the following statement:
With r := sin(x)+x*cos(x),  abs(r) <1/2 is in fact false for some values of x in -1..1. Thus is(abs(r))<1/2 ought to return false. But it returns FAIL meaning in this case (probably) that is cannot determine if it is the case ('inability'  in the quote above).

If you want to use is you must include an elif:

g := proc(x) local r; r := sin(x)+x*cos(x); if is(abs(r) < 1/2) then sin(x) elif is(abs(r))>=1/2 then cos(x) else undefined end if; end proc;

I'm not saying that this makes symbolic integration work in this case (it doesn't), just pointing out the logic.

@Ronan Starting from SHM you can do:
 

SHM := dsolve({Eq1, Theta(0) = 0});
diffSHM := diff(SHM, t);
convert(diffSHM,D);
eval(%,t=0);
eval(%,D(Theta)(0)=Vmax);
solve(%,{_C1});
eval(SHM,%);

The two uses of eval can be combined into one by using eval[recurse]:

SHM := dsolve({Eq1, Theta(0) = 0});
diffSHM := diff(SHM, t);
convert(diffSHM,D);
eval[recurse](%,{t=0,D(Theta)(0)=Vmax});
solve(%,{_C1});
eval(SHM,%);

There are several other ways it can be done, but let me stop here.
 

@9009134 Since dsolve/numeric/bvp uses a finite difference method there is no point in writing your own.

@9009134 The link still doesn't work, but I can see the 3 pages presented as images.
I quote from the first one I can see:
"The three second order ordinary differential equations and their boundary conditions are not in a form particularly amenable to numerical solution."
So I'm afraid I shall have to leave it at that.

@Thomas Dean You need to apply some simplification:
 

restart;
PDE2:=diff(u(x,t),t)=diff(u(x,t),x,x)+2*t*diff(u(x,t),x)+u(x,t)*(ln(u(x,t))^2+ln(u(x,t))-1);
ics2:=u(x,0)=exp(cos(x));
bcs2:=u(0,t)=exp(cos(t^2)),u(1,t)=exp(cos(1+t^2));
sol2:=u(x,t)=exp(cos(x+t^2));
pdetest(sol2,[PDE2,ics2,bcs2]);
simplify(evalc(%));

My guess is that these pdes are made up by basically starting from the solutions. No cheating here because the subject of the paper is not how to solve a pde symbolically.

@9009134 When clicking on the link I get this message:
"The resource you are looking for has been removed, had its name changed, or is temporarily unavailable."
I'm only guessing, but maybe it has been removed by someone for copyright reasons. When I search for the title of the paper I see that it is available for $25.
Anyway, I think that the content of the other worksheet of yours (ich.mw) has been presented by you before. Nothing seems to have changed.

@acer In the help page ?elementwise  we find the statement

"Dimensioned container types: list, set, Array, Matrix, and Vector can be intermixed in a given operation as long as they are the same size."
Right after the example
[true,true,false,false] xor~ <true,false,true,false>;
is given, and another almost identical example with and instead of xor is given in the section Examples.
So in this context, how do you measure size?
Does any version of a list L work on the right hand side of the matrix A here:
 

A:=Matrix(2,symbol=a);
A=~[2,3,4,5]; # error
A=~[[2,3],[4,5]]; # error
V:=Vector(4,symbol=v);
V=~[2,3,4,5]; # Works

## Added:  Pursuing this further (in fact just reading a little more of the help page) we find:
"Lists and sets are always considered to be 1-dimensional. "
I guess that answers my question. But the use of the word size in the earlier statement is too vague to be helpful.