@9009134 I try to keep all my dealings with MaplePrimes matters entirely in this forum. Furthermore, if the reference material you have take up that much space I'm afraid I wouldn't take the time to read it or even just browse it.

@loramina123 You give us h(0)=0 D(h)(0)=tehta (theta, I suppose), D(D(h))(infinity)=1/R and say that infinity might be taken as 2 or 3 micron. By a micron I suppose that (in this context with no units) you mean 1e-6.
But then you also say that you have h(2e-6)=200e-9.  But if infinity can be taken as e.g. 2e-6 then you have a conflict of enormous proportions since R = 1e-7, so 1/R=1e7.
Finally you mention D(D(H)(-infinity)=0 and that -infinity can be equal to 2e-6. I suppose that you meant -2e-6?

In these matters one has to be more careful or else people give up.

@chomchom Are you sure that you want to integrate from -infinity to infinity?

If you are, then try e.g.

display(seq(odeplot(soln1[i], [x, psi(x)], -7 .. 7, color = [red, blue, green][i]), i = 1 .. nops(E)))

If you think that is fine then replace -7..7 by -19 .. 19.

@loramina123 So what are the boundary conditions? You have presented more than one set of conditions. The ode appears to be the same though.
The original 3-point boundary condition problem may be doable, but you do have to use 4 conditions for a 4th order ode. So two conditions could be given at 0, e.g. values of h(0) and D(h)(0). Besides one condition at each of the other two points. By using smoothness at the middle point this can (sometimes) be done.
I did as a rather trivial example the following, where an exact solution can also be found. Thus we can compare the numerical result with the exact.

ode:=diff(y(x),x$4)+4*y(x)=0;
sol:=dsolve({ode,y(0)=0,D(y)(0)=1,y(Pi/2)=0,y(Pi)=0}); #This gives you a simple formula for the solution

By doing what I sketched above I got very good agreement between numerical and exact results.

Before attempting anything similar to your problem I need the actual boundary conditions if you are still interested in the 3-point problem.

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