Bear in mind that I don't really understand what you are trying to achieve!

Your code is a little confusing, because of the way you define ZMAT: it seems to be both a "global" (defined at the top-level) and a "passed-parameter" ( in the procedure definition Z_TABULATE := proc (ZMAT) ).

Since ZMAT only appears to be used within the procedure Z_TABULATE(), I have modified your code so that it is defined and calculated within that procedure. The procedure arguments then become the row/column dimensions. This also means that I don't have to calcuate the row/column dimensions within the Z_TABULATE() procedure

I have made a number of other minor modifications.

The Z_TABULATE() procedure now returns a vector of x-values (PREDS) and the complete matrix ZMAT, which ought to be sets of y-values.

All one then has to do is to create a sequence of pointplot() of PREDS against successive rows (you might want columns??) of the ZMAT matrix

See the attached code: it may not be exactly what you want but It should give you some ideas!

Z_DAK_PROC_Corr.mw

@snhaseela2 

Think of the "mesh" as the 2D grid of points at which the solution to your ODE is calculated.

The exact location of these mesh points can sometimes have an impact on whether a solution can be achieved. One of the most common problems is that if some of the grid points are located along the boundaries of the problem (where your BCs) are defined, there may be enough of a "discontinuity" between the meshpoints along these edges and the adjacent "row" of meshpoints to prevent the numerical solution converging. When this happens you get the error message about the solution being singular at an endpoint. You can (sometimes) get around this error by "offsetting" the grid at which solution points are calculated: this will still use the value along the edges (the BCs), but the actual calculated grid point will be hal-a-step away from the edges. This is achieved using a midpoint, rather than an endpoint method: hence the error message which you received.

As the ?dsolve/numeric/BVP help page says

There are two major considerations when choosing a submethod for a BVP. The trapezoid submethods are generally more efficient for typical problems, but the midpoint submethods are capable of handling harmless end-point singularities that the trapezoid submethods cannot.

The second problem which I had was that the default accuracy (generally set by values of abserr and relerr, which default to 1e-06 and 1e-05 respectively, IIRC) could not be obtained with the default value of maxmesh, which is the total number of grid points being used by the solution method. All (I think?) the solution methods for BVPs use adaptive step control to determine stepsize, but they will not use more than 'maxmesh' gridpoints. Default for maxmesh is 128. Depending on the BVP being solved, the adaptive stepsize control might mean that the required error criteria can be achieved using fewer than 'maxmesh' points. But Maple will not use more. So if you receive the error "unable to achieve default accuracy (abserr=1.0e-06) with the default mesh", you basically have two choices: either slacken the error tolerance (ie increase abserr, relerr) or increase the number of grid points, by increasing the value of maxmesh.

I suspect that the reason Maple does not increase maxmesh automatically, is that you would still need to set some kind of "stopping" criterion in order to prevent the method running for ever on ever-increasing number of grid points. Just think of the current defaults as a compromise, which just "work" for the majority of BVPs.

Since PrebenAlsholm has picked up this problem, my last piece of advice would be to do whatever he advises, because he knows a lot more about solving ODES/PDES than I do!

The sequence of commands

expr:= (a*x+b)/(c*x*x+d*x+e)=k;  # linear over quadratic is equal to a constant
sols:=solve( expr, x);
simplify(subs(x=sols[1], expr));
simplify(subs(x=sols[2], expr));

seems to work OK

The plaintext which you have included is not executable, because it comprises a mixture of Maple input and output. No-one can use this without a *lot* of work, which isn't going to happen. I suggest you upload your Maple worksheet using the big green up-arrow at the right-hand end of the second row in the Mapleprimes toolbars. With this info, we may be able to work on your problem

The sensible interpretation of 1-10 would be -9: now I assumed that you did not mean -9, so I *guessed* what your ambiguous syntax might mean and came up with 1e-10. Just think how much easier life for everyone would be if you used the simple, unambiguous, Maple syntax for a range, as in 1..10

If you look at the output of the code I posted last time you will see that

1. with z=1.0e-10, the solution for Y(eta, 1.0e-10) is a slightly damped oscillation, whose amplitudes start in the range 0.999995..1.000005
2. with z=1.0e-03, the solution for Y(eta, 1.0e-3) is a slightly damped oscillation, whose amplitudes start in the range -50..35
3. with z=1.0 the solution for Y(eta, 1.0) is a slightly damped oscillation, whose amplitudes start in the range -3.3e-53..5.3e53

Now you seem to believe that in somes sense, some or all of these are incorrect - why?

Your second post says

"edit: even though I put z=1e-3 in the initial example it should be more of the order 1-10 and thats where it becomes too large to be believable..."

Your latest post says

"if you set z=1"

so what value of z are you actually interested in? 1e-3, 1e-10, 1 - which part of the word clarify do you not understand?

The code

restart;
vz := 2*(-eta^2+1):
D_im := .22;
r0 := 1:
pde := diff(vz*Y(eta, z), z)-D_im*((diff(eta*(diff(Y(eta, z), eta)), eta))/eta+diff(Y(eta, z), `$`(z, 2)))/r0 = 0:
pde := expand(%):
ibc := [Y(1, z) = 0, (D[1](Y))(0, z) = 0, Y(eta, 0) = 1, (D[2](Y))(eta, 0) = 0]:
sol := pdsolve(pde, ibc, numeric, time = z, range = 0 .. 1):
pds := sol:-value(z = 0, output = listprocedure):
sol:-plot(z =1.0e-10, numpoints = 500, color = blue);
sol:-plot(z =1.0e-3, numpoints = 500, color = blue);
sol:-plot(z =1.0, numpoints = 500, color = blue);

will produce graphs for each of the three cases z= 1e-3, 1e-10, 1. Which of these do you believe to be incorrect, and why?

1. I have no idea where you are getting numbers of 10^50! When I said that the solution was oscillatory - it is oscillating in the range -4..4 (more or less)
2. You say that you expect the solution to be "1 at z=0 and then slightly decreasing as u move along the z-axis", however the command sol:-plot(z = 0.1e-3, numpoints = 500, color = blue); will plot Y(eta,z) versus eta for a fixed value of z=0.1e-3. In other words your are plotting Y(eta, 0.1e-03) versus eta. In this plot, there is no z-axis
3. If I fix eta and plot Y(eta, z) as a function of z, then for all values of eta<1, the solution heads for +/-infinity: with eta=1, then the solution would seem to be well-behaved

Using

eq4 := dsolve({eq2, incs}, y(t), type = numeric, maxfun=0);

the command

plots:-odeplot(eq4, [t, y(t)], 0 .. 5);

executes in 0.327secs (Maple18), 0.390secs (Maple 2015), and 0.436secs (Maple 2016).

On the other hand if I use

eq4 := dsolve({eq2, incs}, y(t), type = numeric, method = lsode[backfull], maxfun = 0):

then the command

plots:-odeplot(eq4, [t, y(t)], 0 .. 5);

takes 10.030 secs (Maple18), 12.370secs (Maple 2015), and 12.932secs (Maple 2016)

So letting Maple select the 'method' is ~30x faster than specifying 'lsode[backfull]' - at least on the three Maple versions I can easily test. (I don't have Maple 17 loaded on my machine, but I very much doubt that the answer would be any different!

My comment about about stepsize control, was based on the fact that when I ran

eq4 := dsolve({eq2, incs}, y(t), type = numeric, maxfun=0);
plots:-odeplot(eq4, [t, y(t)], 0 .. 20);

(which takes about 1.6secs on Maple 2016), I started to see unpleasant "artefacts" on the output plot

These "artefacts" were improved by using

eq4 := dsolve({eq2, incs}, y(t), type = numeric, maxfun=0);
plots:-odeplot(eq4, [t, y(t)],  0 .. 20, numpoints=2000);

(taking 1.65secs) and even further by using

eq4 := dsolve({eq2, incs}, y(t), type = numeric, stiff='true', maxfun=0);
plots:-odeplot(eq4, [t, y(t)],  0 .. 20, numpoints=2000);

which actually reduced the execution time to 0.967secs. The stiff option changes the default solution method from rkf45 to rosenbrock.

eq4 := dsolve({eq2, incs}, y(t), type = numeric, stiff=true, maxfun=0);
plots:-odeplot(eq4, [t, y(t)],  0 .. 100, numpoints=10000);

executed in 4.96secs, and still looked reasonably OK.

Depending on the range over which you wish to plot, and how "pretty" you want the output graph to look, you may be able to improve things further by using any/all of the options maxstep=, abserr= and relerr= in the dsolve() command. From this point, I'm afraid that finding an optimum trade-off between execution speed and acceptable output is a bit of an experimental science

Consider

restart;
f:= unapply
      ( rsolve
        ( { f(n) = f(n-1)+f(n-2),
             f(0) = 0,
             f(1) = 1
          },
          f(n)
        ),
        n
      );
sum( f(n)/n!, n=0..infinity)

@ZAIN UL ABADIN ZAFAR 

The code I sent is syntactically correct.

I have no idea whether it is logically correct - simply because I have no idea what calculation you are performing or what answers you expect

The images which you supplied in your last email do not help much, because I can only actually read one of them

@abdelkarim 

There are still way too many syntax errors for me to work outt what you are trying to achieve. (This is probably no helped by your decision to use 2D math input - since the worksheet you uploaded does not even parse correctly, let alone execute)

I got rid of the parse errors by converting to 1-D input, and fixing the obvious syntax errors - although in many places I had to *guess* what yu were trying to achieve, so some of my "fixes" may well be wrong.

However towards the end of the relevant section, there are some errors(?) in you logic/algorithm becuase I cannot believe that you want what is actuall written - unfortunately I am unable to guess what you do want

The attached worksheet at least executes and you can examine what is being output. I suggest you read the embedded comments very carefully

test.mw

@abdelkarim 

You can upload code using the big green up-arrow in the toolbar for the Mpleprimes editor: right hand end of the second row

There are too many syntax errors in your function definition for me to work out exactly what you are trying to achieve, so I am just going to list a few things which you will have to address

1. Given the starting point (-4) in your 'for' loop, you are trying to evaluate the GAMMA() function at some negative integer values. The GAMMA() function is not defined for negative integers - so this will never work.
2. There is nothing iin your code which is "recursive" - you seem to be generating a straightforward sequence of values, none of which depend on any previous value
3. I'm pretty sure that you do not want to use quantities such as `h[m,t]`. In maple [] are used to denote indexed quantities. So far as I can tell you want 'm' as an index value and t as a simple variable

A quick check suggest that no analytic solution is possible. In order to obtain a numeric solution I calculate that at least four boundary/initial conditions will be necessary as well as numeric values for the indeterminates SI and Pr. Without these, no progress is possible

@adrian5213 

Try the following

randBin2.mw