@dharr 

I almost missed your reply because the contributor counter was stuck at 1.
That is a good tip indeed, thanks

@dharr 

C[1] <= 0: yes of course, what was I thinking?
Only later did I think to check whether the matrix was symmetric.
Finally, I don't understand either why my approach doesn't work after correcting the first condition.
Thanks for you comments

Even if I can't execute your code with Maple 2015, I carefully read it: you used a very beautiful approach and I vote up.

@Rouben Rostamian  

For incormation: I believe this question has connections to this one 
https://www.mapleprimes.com/questions/236730-Solving-Schrodinger-Equation

@somestudent 

It's very tiresome that you spend your time modifying your initial question according to my answers and that you only divulge information on an ad hoc basis.
I feel like I'm wasting my time.

So this is my last answer.
Notice I rewrote your code in a much clever way, for doing copy-paste to define pole1, pole23re and pole23im is a bit stupid (no offense intended).

Your code ends with this command

maximize(MyFilterABS, f = 0 .. .5, location = true)

which is a nonsense for MyFilterABS depends on x, y, w, and f and maximize requires that all the quantities but those you want to maximize with respect to, are given a numerical value.
So I've taken the initiative to replace MyFilterABS by its evaluation when x, y, w given numerical values.

Last point: your sentence "both points are inside the unit circle, both have a maximum at f = 0.5 but for the second one there is one maximum more" is unintelligible.

Download maximum-not-found_mmcdara_2.mw

To find the maxima in a larger f range:

# Ten maxima found

maximize(test_1, f=0..5, location);
numelems(%[2]);
       16.32115142, {[{f = 0.2935930416}, 16.32115142], 

         [{f = 0.7064069584}, 16.32115142], 

         [{f = 1.293593042}, 16.32115142], 

         [{f = 1.706406958}, 16.32115142], 

         [{f = 2.293593042}, 16.32115142], 

         [{f = 2.706406958}, 16.32115142], 

         [{f = 3.293593042}, 16.32115142], 

         [{f = 3.706406958}, 16.32115142], 

         [{f = 4.293593042}, 16.32115142], 

         [{f = 4.706406958}, 16.32115142]}
                               10
# Five maximum found

maximize(test_2, f=0..5, location);
    6.060606061, {[{f = _Z29 + 0.5000000000}, 6.060606061]}
 indets([%]) minus{f}; about(%);
                     {_Z29}
 about(%);
  is used in the following assumed objects
  [_Z29] assumed AndProp(integer,RealRange(0,4))

@apronya 

So much the better

@somestudent 

complexfunction is a function of x, y, f.
Thus the maximizer is a triple (x, y, f=1) and x --> -oo, y --> +oo (obvious result).
The value of complexfunction at this maximizer is equal to +oo (also obvious).

When you write the constraint 

{Maximize(complexfunction, {f <= 1})[1] = 1}

for the "outer" Maximize, you write in fact the constraint

{+oo <= 1}

This is always false and the constraint will never be reached in the "outer" Maximize ... giving you no solution.

I advise you to reformulate your problem correctly before considering solving it with Maple.

Could_it_be_this.mw

@SHIVAS 

On what base do you claim that  what @tomleslie presents is not based on a Finite Difference Method?
It is not because you do not see the code behind pdsolve that pdsolve doesn't use FDM!

Type this in some wotksheet

help(pdsolve[numeric])

In section Details section, you will read this
The default method uses a second order (in space and time) centered, implicit finite difference scheme to obtain the computed solution.

Now look at this page

help(pdsolve[numeric][method])

In section Description it's written :
In addition to the default interface to pdsolve/numeric, which uses the default theta methods (the theta family of methods), additional functionality is available for use of numerical pdsolve for educational purposes.
 ... additional functionality is available for use of numerical pdsolve for educational purposes. [the underlining is mine]
This additional functionality includes the use of 11 standard methods, specification of numerical boundary conditions, specification of startup schemes for two-stage methods, and extensibility for your custom methods

The theta-method is a family of numerical scheme controlled by a parameter "theta".
For 1D-space heta equation the choice theta=1/2 corresponds to the Crank-Nicolson scheme/
Three easy to read references:

• https://folk.ntnu.no/leifh/teaching/tkt4140/._main065.html
• https://www.math.pku.edu.cn/teachers/lizp/courses/Numerical_PDE/Numerical%20PDE%20Lecture%20Slides/numpde_lecture_5_c2.pdf
• https://itp.uni-frankfurt.de/~rezzolla/lecture_notes/2007/ECT_hyperbolic_PDEs_0807.pdf
   

In @tomleslie's code you can force the Crank-Nicolson scheme by setting

 PrVals:=[0.71, 1.00, 3.00, 7.00]:
  colors:=[red, green, blue, black]:
  for j from 1 to numelems(PrVals) do
      pars1:=`union`( pars, {Pr=PrVals[j]}):
      pdSol:= pdsolve( eval([pdes], pars1),
                       eval([conds], pars1),
                       numeric,
                       method=CrankNicholson,
                     );
      plt[j]:=pdSol:-plot( u(y,t), t=2, y=0..inf, numpoints=200, color=colors[j]);
  od:
  plots:-display( [seq(plt[j], j=1..numelems(PrVals))]);
Error, (in pdsolve/numeric/par_hyp) the CrankNicholson method can only be used for a single PDE

The error message is clear.
I guess (see above) that the development team didn't bother to implement the Crank-Nicolson scheme for pde systems for two reasons:

• Crank-Nicolson scheme is mainly for educational purposes.
• and theta-method is more general.
  So why do the job twice?

 Note that if you really want to implement the Crank-Nicolson scheme for pde systems, you can do it yourself, see here

help(pdsolve[numeric][pdemethods])

As your question was "How to get same graph from maple with finite difference method for differential equations " I consider that @tomleslie has perfectly answered it (and I vote up).

Could you provide us the expressions of g and f?
It is absolutely necessary for one can build extremely simple examples where no solution does exist:

Download SolutionNotGuaranted.mw

@AHSAN 

You can customize "my" code the way you want to get the rendering which pleases you.
But writting a procedure which accounts for all your desires is, IMO, a total waste of time.

Nevertheless here is an evolution of my previous code.
 

Download Evolution_1.mw

 

If tou wnat to change the indication at he end of a bar, you must add an optional parameter to procedure UnivariateBarGraph (all these optional parameter are gathered in curly braces in the sequence of parameters of the procedure).
For instance:

 , HeadTitle::string:=""

and change

plots:-textplot(
        [
          values[i]+Hpos
          , nuval[i]
          , nprintf("%1.3e", nuval[i])
        ]
        , align=`if`(
                   Halign = NULL
                   , `if`(nuval[i] > 0, above, below)
                   , `if`(nuval[i] > 0, {above, Halign}, {below, Halign})
                 )
      )

by

plots:-textplot(
        [
          values[i]+Hpos
          , nuval[i]
          , cat(HeaderTitle, sprintf("%1.3e", nuval[i])
        ]
        , align=`if`(
                   Halign = NULL
                   , `if`(nuval[i] > 0, above, below)
                   , `if`(nuval[i] > 0, {above, Halign}, {below, Halign})
                 )
      )

@AHSAN 

... the most dramatic is that you must not assign beta to a value or a list where tou did that.

Look to the attached file
 

Download UnivariateBarGraph.mw

@AHSAN 

Nu is a function of Q, beta and lambda:

indets(Nu, name);
                       {Q, beta, lambda}

You write: "I want to draw bar chat for different value of beta [0.1,0.3,0.5,0.7,0.9) between Nu and y ( Nu is on vertical axies and y is on horizontal axis) and crosspedning to the values of beta the value of Lambda and Q also changes".

What is y (which is not part of the definition of Nu) ?

One doesn't draw (a) bar chart for different value of beta  between Nu and y... one draws (a) bar chart of N (and y if any) for different value of beta : I still do not understand what you want.

It would be helpfull if you coukld provide a sketch of what you want (use PowerPoint or do this drawing by hand and scan it).

The content of your mw file has nothing to do with data that can be presented under BarChart form.
Basically it contains an expression N which depends upon beta, lambda, Q and that's all.
See here 

help(BarChart);

Maybe you should explain clearly what your desired BarChart is aimed to represent?

@arashghgood 

Could you provide your Maple code?
Just click on the big green up-arrow in the menu bar. 

@acer 

Thanks acer for all these precisions, particularly for the trick on your point (2).

@Carl Love 

Well seen.
I always try to be very careful with these Hfloats, which are automatically switched on with some (all) of the functions in the Statistics package (the function to plot is a scaled PDF of a ChiSquare(N)^2 random variable with large N).

UseHardwareFloats:= false:  works perfectly well.

Thanks