@tomleslie I tried making a new document in document mode and 2d input, just 3 lines. I copied this and pasted it into a fresh worksheet (worksheet mode and 1d input).
There I converted the pasted code into 1d-input. That worked, but there are no prompts or execution group delimiters and the area covered by the lines seems to be a document block.
This can be changed by using the menu item Format/Remove Document Block.
Unfortunately this meant that (in my case) 3 lines of
print(??); # input placeholder
also appeared.
#############
I just had a look at the discussion about converting a file made in document mode and 2d input into a file in worksheet mode and 1d input.
http://www.mapleprimes.com/questions/203080-How-To-Convert-Document-In-Document

I had occasion to do this the other day
http://www.mapleprimes.com/questions/204487-Too-Lengthy-Parametric-Equations

It would have been easier to do the conversion by exporting to .mpl as suggested by Alejandro Jakubi in the first reference.

@cuongtd The output from allvalues(solve....) is not a numerical result, but is the exact result.
As it turns out all 32 solutions (16 real and 16 imaginary) can be expressed in terms of the 16 roots of a polynomial of degree 16. Once you use evalf on that exact result numerical rootfinding is used, but since the problem just lies with the polynomial surely all roots are found.
To see the result before using evalf just change the colon in the first line below into a semicolon: It is huge!
#I'm assuming that eqs and eqs2 are defined as I did in my answer. Still using Digits=30.

resA:=allvalues([solve(eqs union eqs2)]): ##Sequence of length 16 each consisting of a list of two solutions
nops(resA[1]); # A pair of 2 solutions
nops([resA]); #16 pairs
nops(ListTools:-Flatten([resA])); #32 solutions
indets(resA[1],specfunc(RootOf)); #Notice the two arguments to RootOf: a polynomial of degree 16 and an index
ROF:=indets([resA],specfunc(RootOf)); #All of the RootOf's
nops(ROF); #16
pol:=map2(op,1,ROF); #The same polynomial for all
map2(op,2,ROF); #The 16 indices
#Now using evalf as I originally did.
res:=evalf(resA):
resR:=remove(has,[res],I); #Looking at the real results only
##Comment: It so happens that each two in the 16 pairs are either both real or both imaginary. I checked that.
map(nops,resR); #8 pairs
resRflat:=ListTools:-Flatten(resR);
nops(%); #16
##Graphical illustration: Animation showing one solution at a time.
##The 4 points in each picture are (ci,si), i=1..4.
plts:=seq(plots:-pointplot(subs(resRflat[j],[seq([cat(c,i),cat(s,i)],i=1..4)]),color=[red,blue,green,maroon],symbolsize=25),j=1..nops(resRflat)):
plots:-display(plts,insequence);
##Note: Recall that my (c4,s4) is the point (cos(t4),sin(t4)), where t4 = t+t1.

@mapleus I notice that A(0) is imaginary. Is that intended?

eval(A(0),r(0)=0.3-0.01);
#Outout: 9.330399019*I

Before trying a loop like that it makes sense to make it work for just one step!
That is, make dsolve come out with some result. I got an error.

Besides that I see one problem right away:
subs(F,r(s)) is itself a procedure, so you should just have
r_ravn:=subs(F,r(s));

@nafiqah38 As you can see yourself there is nothing at all in your reply.

Give us your system with boundary conditions etc. either as pasted code (text) or as an uploaded worksheet. For the latter use the thick green arrow in the editor.

@FMEHR The approach I showed above works with an appropriate extra condition.

restart;
dsys3 := {-901.7765286*(diff(f1(x), x, x))+3.153274902*10^(-10)*(diff(f2(x), x, x, x))+0.4807756008e-5*(diff(f2(x), x))-27.92641797*(diff(f3(x), x))+18074.25056*f1(x)-3.100642108*(diff(f3(x), x, x))*(diff(f3(x), x))+7.74301535*(diff(f3(x), x))*f3(x) = 0, 0.1873004296e-2*(diff(f3(x), x, x, x, x))+(diff(f3(x), x, x))*(-0.1142850079e-2-0.1587291776e-1*P)-.1797931318*(diff(f3(x), x, x))+13.96314014*(diff(f1(x), x))+1.822500000*10^(-10)*(diff(f2(x), x, x))-4.888800000*10^(-7)*f2(x)+22.72337812*f3(x) = 0, 0.2249991345e-1*(diff(f2(x), x, x, x, x))-227.1861332*(diff(f2(x), x, x))+29537.78948*f2(x)-2.841067252*10^(-11)*(diff(f1(x), x, x, x))-3.39267254*10^(-7)*(diff(f1(x), x))+7.650000000*10^(-11)*(diff(f3(x), x, x))-1.400000000*10^(-7)*f3(x)-5.500000000*10^(-7)*f3(x)^2+9.278333333*10^(-13)*(diff(f3(x), x))^2+3.300000000*10^(-13)*(diff(f3(x), x, x))*f3(x) = 0, f1(0) = 0, f1(1) = 0, f2(0) = 0, f2(1) = 0, f3(0) = 0, f3(1) = 0, (D(f2))(0) = 0, (D(f2))(1) = 0, (D(f3))(0) = 0, (D(f3))(1) = 0};
 
indets(dsys3,name);
for i to 3 do dsys3[i] end do;
diff(dsys3[2],x);
eq:=subs(solve(dsys3[3],{diff(f2(x),x$4)}),%);
sys:=solve({dsys3[1],dsys3[3],eq},{diff(f1(x),x$3),diff(f2(x),x$4),diff(f3(x),x$4)});
dsys3[4..-1];
cdn:=eval(convert(dsys3[2],D),x=0);
#Using the extra condition f3'''(0)=-1 we have success with
res:=dsolve(sys union dsys3[4..-1] union {cdn, (D@@3)(f3)(0)=-1},numeric);
res(0); # Eigenvalue P = -42.85...
#If we plot the three together two of them are drowned by the largest (f3). So scale:
plots:-odeplot(res,[[x,100*f1(x)],[x,10^11*f2(x)],[x,f3(x)]],0..1,thickness=3,legend=[400*f1, "10^11*f2", f3],size=[1024,600]);
#Or in 3 separate graphs
use plots in
  display(Array([seq(odeplot(res,[x,cat(f,i)(x)],0..1),i=1..3)]))
end use;
## Finding other P's
## You can try changing 'a' in the extra condition (D@@3)(f3)(0)= a (a<>0). I tried a = 0.4 as an example (and started with Digits:=30). In the linear case this change should not be expected to change matters a priori since if f = (f1,f2,f3) is an eigenvector (function) so is a scalar multiple of f. But your system is nonlinear.
I got P = -41.10... and the graphs were different.

Take a look at the thread:

http://mapleprimes.com/questions/204463-Nonlinear-Ode-Equations

where someone else asked a very similar question just using different constants, but the very same notation: dsys3, f1,f2,f3.
My guess is that it is a fellow student. Talk to him.

@emmantop What I meant is illustrated in this very simple example:
restart;
ode:=diff(x(t),t,t,t)+x(t)/(3-t)=0;
dsolve({ode,x(0)=0,D(x)(0)=1,x(3)=0},numeric); #Error because of the denominator 3-t.
res:=dsolve({ode,x(0)=0,D(x)(0)=1,x(3)=0},numeric,method=bvp[midrich]);
plots:-odeplot(res,[t,x(t)],0..3);
res(3); #Here we get the values of x, x', x'' at 3
res(0); #Here we get the values of x, x', x'' at 0

@FMEHR The image above was produced by my own unsophisticated solver. I have not yet been able to get that by dsolve.

After observing that it produced a result with f2 = 0 I tried using that directly as follows (assuming that sys and cdn are defined as above):

sys20:=eval(sys,f2(x)=0);
for i to 3 do sys20[i] end do; #Viewing the result
cdn20:=eval(cdn,f2=0);
bcs20:=remove(has,dsys3[4..-1],f2);
#Leaving out sys20[1]:
res0:=dsolve(sys20[2..3] union bcs20 union {cdn20},numeric,approxsoln=[f1(x)=-sin(Pi*x),f3(x)=100*sin(2*Pi*x)]);
res0(0.5);
plots:-odeplot(res0,[[x,f1(x)],[x,f3(x)]],0..1); #Pretty small: Reliable???
#Checking if sys20[1] happens to be satisfied:
plots:-odeplot(res0,[x,lhs(sys20[1])],0..1,thickness=5); #It seems so

@FMEHR Do you have any reason to believe that a nontrivial solution exists?

Possibly something like the following, where f2 is the blue graph, i.e. f2(x) = 0 (almost at least)?

I don't trust the solution above since I cannot make dsolve come out with that answer even if I feed it as an approximate solution. It still returns the trivial solution.

All your constants need to have numerical values. So we need to know what those are.

The boundary conditions you give as:

bcs := y[1](0) = 0, y[2](0) = 1, y[3](0) = d[1], y[2](blt) = 0, y[3](blt) = 0, y[5](0) = 1, y[5](blt) = 0, y[4](0) = d[2], y[6](0) = d[3], y[4](blt) = d[4];

where blt is known, but d[1],..., d[4] are not. The six boundary conditions not containing the d's should determine the solution, so the values of y[3](0), y[4](0), y[6](0), y[4](blt) will be determined directly from the solution.
This has nothing to do with using a shooting method or not.

@emmantop Surely, you had a system of two odes of order 4 in f and order 2 in theta.
For some reason (which?) you want to turn this into a first order system using the variable names y[1],...,y[6].
Thus you must write down 6 first order equations of the form
diff(y[1](eta),eta)=y[2](eta);
diff(y[2](eta),eta)= whatever it is in terms of all the y[i]'s.
etc. up to and including
diff(y[6](eta),eta)= whatever it is in terms of all the y[i]'s

You have only two of these at the moment, but you have a line containing the translations from f,theta to the y[i]'s.
So the information from there has to be written in the form I showed above.

Basically, DEtools[convertsys] could have done that.
Still the question remains, why convert to a first order system? dsolve/numeric will do it itself (without telling you).

A suggestion if you insist on writing the original system of 2 equations (say eq1, eq2) in f and theta as a first order system:
Write down eq1 and eq2. Then do
DEtools[convertsys]({eq1,eq2},[],[f,theta],eta,Y,YP);
Then you get a list, whose first element is the first order system. Now write that in the form desribed above.

Before getting to the issue at hand you need to fix (at least) two things:

1. Missing multiplication sign in ode1 after
2. There is a theta(eta) appearing in ode2. That should probably be y[5](eta).