@loramina123 As Tom Leslie has examined it is difficult if not impossible to solve the revised problem by dsolve/numeric/bvp.

The original problem involved conditions at 3 points, say a,b,c. Was that the real problem you wanted to solve?
Certainly it cannot be done in a straightforward way by dsolve/numeric/bvp, but maybe you could split the problem into two boundary value problems for the same ode but on the two intervals a..b and on b..c.
Extra conditions at b would be the requirement that the solution be smooth (also) at b.


Do your revised conditions also reflect a change in the actual (physical) problem?

Your code is extremely difficult to read since it is probably some copy of 2D-math input.

But you cannot use psi(infinity)= whatever as a boundary condition in dsolve/numeric/bvp. You must replace infinity by a real number.

@YasH I tried the following variant of your code:

for k from 1 to 4 do [Transpose(Eigenvalues(M)),Transpose(op(1, [Eigenvectors(M)] ))] end do;

and found that the eigenvalues in the first version of the eigenvalues vector (the one from Eigenvalues) always came in the same order for the given matrix M whereas this was certainly not the case for the second version.
I have no idea why.

@9009134 Do you have a reference for this problem? I should like to see how SYS was derived.
I notice that the coefficients in SYS are either 1, -1, 1/2, or 3/2, i.e. there aren't any values that come from physical data as opposed to quite a lot of bvp problems presented on MaplePrimes.

It is safe to say that none of us like images if we have to examine some problem. We want code in text or an uploaded worksheet: Use the big fat green arrow in the MaplePrimes editor to upload a worksheet.

@9009134 My point in trying revised boundary conditions was actually to make you examine the physical situation of which (presumably) this ode bvp is a mathematical model. Could the boundary conditions be different; why impose f''(1)=0? Is there a good physical reason?
Obviously, I wouldn't have tried revising the bcs if we hadn't run into problems like hints at nonuniqueness with the original bcs.

For the fun of it (seriously, yes) you could try deactivating the error check by setting abserr ridiculously high.
It won't hurt to try. Here I'm using abserr=10^10, which basically means that that all components of the corresponding first order system are allowed to stray 10^10 away from the "true" solution!!!
You won't be allowed to use abserr = infinity, but 10^10 is just a substitute.

Digits:=15:
res3 := dsolve(sys2 union bcs3, numeric, maxmesh = 2048, abserr = 10^10,initmesh=256);
plots:-odeplot(res3, [[y,g1(y)],[y,100*g2(y)]], 0 .. 1);


Whether this has anything whatsoever to do with a solution I don't know.
As a warning you may try plotting the derivatives of g2. Here the third (corrected from 2nd) derivative:
plots:-odeplot(res300, [y,diff(g2(y),y,y)], 0 .. 1);

@9009134 I forgot to give the plotting command that produces the graph shown:

plots:-odeplot(solA,[seq([eta,fu(eta)],fu=[10^3*F,10^6*K,10^8*Omega,theta])],0..1,legend=[10^3*F,10^6*K,10^8*Omega,theta],color=[red,blue,green,black]);

As I said I got inspiration to the approximate solution from results obtained with my private program. These results were in matrix form. They were used as an approxsoln in dsolve/numeric/bvp and the graphs in my comment (the first set) were produced. Then second time around I simply looked at those graphs and took out of my head the approxsoln, which basically was
[F(eta)=eta^2,K(eta)=tanh(eta),Omega(eta)=tanh(eta),theta(eta)=1-eta]
but subsequently modified with scaling factors to make the graphs of those look like the graphs in my comment.
Thus I wasn't doing anything fancy this second time.

It is disturbing that the second set of graphs (the graphs in the "answer") although having the same form are clearly different from the set of graphs in my "comment", e.g. does the graph of 10^6*K level off at about 2.5 in the "answer", but at about 1 in the "comment". Nonuniqueness is possible with your boundary conditions, but I wouldn't bet much on the results.
Please see,however, the addition to my answer where I replace the requirement (D@@2)(F)(1)=0 by F(1) = 1.
##
In the help page for dsolve/numeric/bvp it also clearly states that the methods used are not well suited for problems with singularities.

@Naeem Ullah I don't see a link to the worksheet in your latest reply.

@Naeem Ullah But you need to give that gamma2 a value. What would it be? The value will obviously influence the results.

By a greater value for abserr I simply meant a value greater than 1e-6, as e.g.  1e-3 .
You may even go back to "inactivating" the error check by setting it ridiculously high like you had originally (10^20) just to see what you get.

It is worth mentioning that the help page for dsolve/numeric/bvp says this about the methods used:
... the midpoint submethods are capable of handling harmless end-point singularities ...

and continues with

The available methods are fairly general, and should work on a variety of BVPs. They are not well suited to solve BVPs that are stiff or have solutions with singularities in their higher order derivatives.

As for contact through email I prefer to discuss the issues in this forum.


 

@Naeem Ullah You said that gamma is the curvature parameter and gamma1 is just a constant set at 1.

But gamma in your new version of the worksheet is still just Euler's constant, and is an initially known constant in Maple.
Thus it couldn't be a curvature parameter as it appears now. It will be treated as a constant with value approximately  .5772156649 as I said earlier.

You never assign a value to gamma, and in fact you won't even be allowed to:
gamma:=1.2345; # results in the error:

Error, attempting to assign to `gamma` which is protected.  Try declaring `local gamma`; see ?protect for details.

As the error message says: You can start your worksheet like this:
restart;
local gamma;
##Then the rest. Now you can assign to gamma.
I find it simpler just to use another name instead of gamma, e.g. gamma2.

About abserr: I didn't mean to say that you have to use the default abserr, which is 1e-6. You may have to raise it.
I was just pointing out that abserr=10^20 seems ridiculous as it would mean that all that is guaranteed about the solution is that it hopefully would not be more than 10^20 away from the true solution!
It does appear, however, that dsolve/numeric/bvp runs through some initial work which may bring forth a solution that is actually not bad, but you will have no guarantee that it is any good.
We may conclude that using abserr=10^20 or some huge number like that basically makes dsolve/numeric/bvp skip error checking.
Skipping error checking isn't all that bad an idea if you are desperately trying to get some idea of how the solution looks.

@Naeem Ullah You have Euler's constant gamma ( evalf(gamma) = .5772156649) in some places and gamma1 in other, where gamma1 is later defined to be 1. Is this intentional or should both be gamma1?

It is rather strange to me to see abserr=1e20 considering that the default is 1e-6.

Clearly you don't use or need the deprecated linalg package. If you need LinearAlgebra somewhere later, use LinearAlgebra.

Do you really want to set Digits:=24?

@nm As one user who has answered (or tried to answer) an awful lot of questions about boundary value problems for third or fourth order odes I must agree with Carl; there have been very many, but most have had to do with convergence problems, not with the syntax for higher derivatives in the bcs.
Looking under my own name in Users, All answers, first page, the first question about a bcs problem is this one in which the second derivative happens to appear in the bcs:
http://www.mapleprimes.com/questions/218225-Error-Showing-As-#answer232183

@artfin I wasn't referring to your use of t->...., but to the first argument to Int, which wasn't a procedure, but was a procedure call, i.e. sin(x) versus sin.

The t from the error doesn't come from the obvious t in the outer construction of yours:
comp_1:=t->evalf(Int(x->r1_lab_deriv(x)[1],0..t,epsilon=1e-6,method=_d01akc)):

as will be seen if you change that t to e.g. s (in both places). The error still refers to t.

The error is due to the t in r1_lab_deriv and inside that it is due to the term r1_rot_deriv(t).
r1_rot_deriv(t) you have defined as
r1_rot_deriv:=x-> <fdiff(r1_rot(t)[1],t=x),fdiff(r1_rot(t)[2],t=x),fdiff(r1_rot(t)[3],t=x)>: using fdiff.
Furthermore
r1_rot(t) is defined as
<-r1(t)*cos(q(t)/2),0,-r1(t)*sin(q(t)/2)>:
where r1 and q are procedures returned by dsolve/numeric (result in solution_procedure).
The error comes from the fact that procedures returned by dsolve/numeric with output listprocedure return unevaluated when receiving just a name like t.
Since it follows from eqs that diff(r1(t), t) = 0.5443205700e-3*p1(t) you could define
r1_deriv by 0.5443205700e-3*p1 and similarly q_deriv by  using the odes.
The less convoluted the end procedure is the better. You have sequences of procedures calling each other.
##
But it may be better (and easier) to include needed procedures in the original call to dsolve/numeric.
I shall confine myself to a simple example:
##
restart;
ode:=diff(x(t),t)=sin(x(t))+t;
res:=dsolve({ode,x(0)=0},numeric,output=listprocedure);
X:=subs(res,x(t));
fdiff(X,[1],[1.2345]);
eval[recurse](rhs(ode),{t=1.2345,x(t)=X(t)});
## Now include the derivative as xp instead:
res2:=dsolve({xp(t)=rhs(ode),lhs(ode)=xp(t),x(0)=0},numeric,output=listprocedure);
X,XP:=op(subs(res2,[x(t),xp(t)]));
XP(1.2345);

@Naeem Ullah There is no way we can help unless we get the entire code, either as text (not images) or as an uoloaded worksheet. Use the fat green arrow in the second line to the right in the MaplePrimes editor.
We need to be able to work with the code ourselves without having to type in your code first.