Clicking on the link returns "page not found". Likely because of non latin characters in its name.

@AHSAN 

        CODE REMOVED BY THE AUTHOR  see here

ALL CODE SNIPPETS DELETED BY THE AUTHOR  see here

[1]
What does diff(p(x), x)  mean in equation 5?
Wouldn't it be diff(p(y), y)  instead? But if it is so you cannot solve equation 5 numerically as p(y) is undefined.
  

[2]
You have 2 odes (equations 5 and 6) that psi(y) must verify... you must decide the one you want to keep unless they both have the same solution, whixh is highly unlikely given thefirst one depends
on p.

[3]
Your bc are given at 3 points (-1, 0, 1) which is not supported by Maple.
Either you formulate the problem differently by considering a 1st problem in [-1..0], a second in [0..+1] and interface continuity equations... or you remove on point, for instance 0


[4]
I psi(y) verifies equation 6 and bcs 7 you can plot u(y) this way

        SNIPPET DELETED BY THE AUTHOR  see here

And if you want to plot u(y) for different values of k:

​​​​​​​[5]I psi(y) verifies equation 6 and bcs 7 you can plot the rhs of de1 this way
  
​​​​​​​

If it is p(y) that you want to plot, its expression is given at the end of point [3] for k=0.


Reformulate your problem correctly and come to us when it's done

@C_R 

I found this inadvertedly and I use it a lot, but you are right, it is likely an undocumented feature.
I guess the option is accepted because its syntax is the same as PLOT's

Compare these:

plot(1, x=0..1, color=blue); # documented
op(%);  #option 1

plot(1, x=0..1, color=COLOR(RGB, 0, 0, 1)); # undocumented?
op(%); #same structure than option 1

plot(1, x=0..1, color=ColorTools:-Color([0, 0, 1]));  # documented in plot/color
op(%); #same structure than option 1

@vv

It seems that the orthogonal polynomials which ease the interation the most are Chebyshev U.
If the solution of order n is S_n(x) = C₀*U(0, x) + ... + C_n*U(n, x), then the solution of order (n+1) should be 
S_n(x) = S_n(x) + C_n+1*U(n+1, x).

More of this the results are significantly better when numeric integration is done instead of exact integration(???)

Galerkin_solution.mw

@vv 

In fact, I had something else in mind that I didn't see through to the end.

Using othogonal polynomials enables building successive approximations but just computing a single new term for each refinement.
So I began coding this with Hermite polynomials in mind, while I realized that the singulatity of

phi(y)/(x-y+2)^2

will cause problems to compute 

Int(H(m, y) / (x-y+2)^2 *  exp(-y^2), y = -infinity..+infinity)

So I resricted th integration range in order to avoid this singularity but at the same time I unfortunately kept the Hermite poynomials.

I might as well have written

Galerkin := proc(EQ, M, R, d)
  local G, S, N, r, m, n, C:
  G := Matrix((M+1)$2):
  S := Vector(M+1):

  r := evalf(R):

  Digits := d:
  for m from 0 to M do
    N := int((x^m)^2, x=r);
    for n from 0 to M do
      G[m+1, n+1] := int( (value(eval(lhs(EQ), phi = (u -> u^n))) assuming x > 0) * x^m, x=r) / N
    end do:
    S[m+1] := int(rhs(EQ)*x^m, x=r) / N;
  end do:
  G := evalf(G):
  S := evalf(S):
  C := LinearAlgebra:-LinearSolve(G, S):

  return add(C[i+1]*x^i, i=0..M)
end proc:

which gives (obviously) the same final result (R=0..1 for instance)

Thank for pointing out this error.

@Carl Love 

but it's too top-notch programming for me to be able to appreciate all its subtleties.

@Carl Love 

You use gamma in the third edo but gamma is never defined.
By default gamma is the Euler constant (~0.5772156649).
So I advice you to give gamma a correct value (use local gamma at the top of your worksheet) before applying dsolve.

Loot to what happens here

dsolve(eval({ics, couple[]}, gamma = 10^9), numeric, maxfun = -1)

PS: had you use another name than gamma , dsolve would have sent you this type of warning:
Warning, The use of global variables in numerical ODE problems is deprecated, and will be removed in a future release. Use the 'parameters' argument instead (see ?dsolve,numeric,parameters)
This is because gamma is a "predefined" name that you did not get it.

@Carl Love 

However, I find it a pity that orthopoly is present in the list of packages without any warning as to how to use it within a procedure.
Would it really be that complicated to rewrite orthopoly as a module?

Even if your integral equation is a convolution equation, usually handled using Laplace transform, it seems that this approach does not provide a result.

In the attached file I use the Galerkin method to build successive approximate solutions over the range x=0..1.

Galerkin_solution.mw

For instance

phi(x) = .250008-.124777*x+0.933750e-1*x^2-0.708356e-1*x^3+0.493446e-1*x^4-0.187361e-1*x^5+0.100280e-1*x^6-0.250550e-1*x^7+0.155970e-1*x^8

verifies your equation with an absolute error less than 7.5734e-05

Download Of_course_it_is_correct.mw

Given the error message you get I would say you di not use fsolve correctly. Look at this:
 

Download fsolve.mw
 

@Carl Love 

"Unit is exactly the same command as Units:-Unit"
In my defense, I have to say that the help page isn't very helpful in this respect:

Thank you for drawing this to my attention.

@emendes 

I suggest you to identify the function names BEFORE applying ddf (or df):

# Your last question... first way:

PS1 := {
         Phi(z,n=1)*Phi(y,n=1)*Phi(y,n=2)+3*Phi(x,n=2), 
         Phi(z,n=1)*Phi(y,n=2)+6*x*y,x*Phi(y,n=1)-2*Phi(x,n=2)*y
       }:

indets(PS1, name) minus {n}; 
              {x, y, z}

eval(PS1, Phi=ddf);  

In my opinion this is the simplest way to proceed.

But you can also "extract" the indeterminates from the output of PS1 (even if it is more complex)
 

# Your last question... second way:

PS1 := {
         ddf(z,n=1)*ddf(y,n=1)*ddf(y,n=2)+3*ddf(x,n=2), 
         ddf(z,n=1)*ddf(y,n=2)+6*x*y,x*ddf(y,n=1)-2*ddf(x,n=2)*y
       }:

i  := indets(%):
si := convert~(i, string):

use StringTools in
  map(ii -> if Search("mi", ii) = 0 then parse(ii) end if, si);
  map(ii -> if Search("mi", ii) <> 0 then parse(substring(StringSplit(ii, "mi")[2], 3..-4)) end if, si);
  print('indets' = % union %%)
end use:
               indets = {x, y, z}

An_alternative.mw