@assma OK, so wr(x,y) and w3(x,y) are functions.  But in your original post you proposed to compute their difference through:

    LL := proc (x, y)
      option inline;
      abs(wr(x, y)-w3(x, y));
    end proc;
    maximize(LL(DD*i, DD*j), i = 0 .. N, j = 0 .. N)));

What is the reason for the i and j grid?  Why not simply do:

    maximize(LL(x,y), x=xmin..xmax, y=ymin..ymax);

Edgardo, I downloaded your latest Physics package.  Tried out the similarity solution and saw that it works as intended.

I sincerely thank you for your contributions to Maple in general, and for your immediate attention to users' input in particular as we have seen repeatedly over the years.

Rouben

@Jawadqau That raises several more questions.

• What do you mean by u=crf?  On the left we have u(r,z).  On the right we have c*r*f(eta).  How is z related to eta?
• From w=(c*nue*h)^1/2  it appears that w(r,z) does not depend on r.  Is that correct?
• Apparently (u,v,w) are velocity components, but what in which coordinate system?  The notation u(r,z) gives the impression that the coordinate system is cylindrical.  Are (u,v,w) components in the cylindrical system or in an associated Cartesian coordinate system?  If cylindrical, then probably u is the radial component of the velocity and v is the tangential component.  But that's only a guess.
• Are you asking for streamlines before or after the similarity transformation?  These will be very different since the transformation w=(c*nue*h)^1/2  is nonlinear.  To plot the streamlines in the original variables, you need to supply numerical values for all the coefficients that enter the transformation.

@assma Let me repeat: Forget about the CPU time for now. Check your code for correct logic.  The little that you have shown does not look correct.

What is w3(x,y)?  You have said it is an approximation.  Can you say a few more words? What sort of approximation?  What is the underlying problem that you are attempting to solve?

Is it a finite difference approximation?  Yes or no?

If it is a finite difference approximation, then it must be defined on a discrete grid, therefore w3 should be an array, not a function.  That's why I suspect that your use of w3(x,y) is fundamentally wrong.  Perhaps you can explain your reasons.

We will talk about the CPU time after you have clarified these issues.

@assma First, let's verify that you have a working code.  We will worry about its CPU time later.  You haven't said whether your code works.  Does it?  Do you get a numerical value for the L-infinity error?

@Jawadqau Did you forget to tell us how f, g, h are related to u, v, w?

@assma OK then, I assume that you are telling me that your question has been adequately answered.

@assma Here is how:

err := abs(wr(i*dx,j*dy) - w3(i*dx,j*dy));
max_err := max([seq(seq(err, i=0..N), j=0..N)]);

What I have shown here is a direct answer to your question. I suspect, however, that you are not formulating your question correctly, therefore this answer may be of no use.

Specifically, it seems to me that you wish to calculate the error in a finite difference approximation of a differential equation.  If so, then it is likely that w3(x,y) is not  a function at all.  Rather, it is an array w3[i,j].   I cannot provide more help with the little information that you have provided.

What are wr() and w3()?  What is N?

@waseem The two calculations you have shown are not comparable; you can't do in one as you do in the other.

The calculations in your first post express lambda (vertical coordinate) as a function of delta2 (horizontal coordinate).

The calculations in your second post express V (horizontal coordinate) as a function of y (vertical coordinate).

Note the switching of the horizontal/vertical.

@Alex Bowden 

The mathematical notation (a+b)(x) is inherently ambiguous—is that a multiplication or function application?

Humans resolve the ambiguity in context as vv has already pointed out.  The context is mostly in one's head.

Absent that context, Maple's 2D input attempts to resolve the ambiguity by applying certain rules.  I don't think that it is possible to design rules that satisfactorily resolve all cases.  Not simple rules, anyway.

Maple's 1D input provides a perfectly context-free solution.  Many people, including me, prefer the 1D input for that reason.  You may want to give it a try.

@mehdibaghaee In the first example, two planes in 3D intersect to create a line.  The calculations lead to a parametric representation of that line.  Eliminating the parameter makes no sense. The parameter is an essential ingredient of the equation of the line. The equation y+2*x=7 that you have suggested is not the equation of a line at all; it is the equation of a plane that contains the line.  Having  eliminated the parameter, you have lost the line.

@ecterrab 

I take exception to the statement: "look at the solution, with the term under focus of the form _C2*erf(abs(x)) where x::real.  You know erf(-x) = -erf(x); reabsorb now that sign within _C2"

There's the rub.  You cannot absorb that sign into _C2 because the derivation of the solution presupposes that _C2 is a constant.  Now you are making it depend on x (in the form of a step function).

To see the error in that argument more clearly, let's look at another elementary example.  Unlike my first example which was an initial value problem of which you did not approve, this one is a boundary value problem, just like that in the original post.

@ecterrab 

I agree with most of what you wrote, but I disagree on the "this is not really a bug" part. A function that satisfies a differential equation everywhere except at one point, is not a solution, and a program that produces such a function is buggy and needs to be fixed.

To make this absolutely clear, let's look at this very elementary example.

@torabi It is not difficult to show that zero is the only possibe answer.  Here is how.