@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).

@surnizai The plot of the complex valued function x*y1 (defined in the complex x-plane,but restricted to the rectangle x=-2-I..2+I) is really a 3d-plot with the "x"-axis being Re(x), the "y"-axis being Im(x), and the "z"-axis being the absolute value of the function, i.e. abs(x*y1).

@Carl Love I tried the following:

res:=dsolve(eq);
ode:=op([2,2,1,1],res);
resb:=dsolve(eval(ode,epsilon=0));
eq3:=t=subs(_a=x,int(1/rhs(resb[1]),_a))+_C2;
eq4:=t=subs(_a=x,int(1/rhs(resb[2]),_a))+_C2;
sol1:=solve(eq3,x);
nops([sol1]);
simplify([sol1]); #Not simple
# It turns out that the other solution for b(a) leads to the same two solutions for x:
sol2:=solve(eq4,x);
nops([sol2]);
sol1[1]-sol2[1];
sol1[2]-sol2[2];

I don't know how to make it as simple as the direct answer from dsolve with epsilon=0.
But introducing two other arbitrary constants replacing _C1 and _C2 might get us somewhere closer.
indets(sol1[1],radical);
#This example is somewhat encouraging:
evalf(eval(sol1[1],{_C1=1,_C2=0}));