restart:
eq := (-2*Psi^2*(a[0]*(Psi-1)^2*(Psi+1)^2*cosh(d)^4-(Psi+1)*(Psi-1)*(Psi^2*a[0]+(-a[1]+b[1])*Psi-a[0])*cosh(d)^2+Psi*(Psi^2*b[1]-Psi*a[0]+a[1]))*abs(a[1](t)*Psi(t)^2-a[0](t)*Psi(t)+b[1](t))^2+((Psi+1)*(Psi-1)*(I*Psi^4*a[1]-2*Psi^3*a[1]+(-I*a[1]-I*b[1])*Psi^2-2*b[1]*Psi+I*b[1])*cosh(d)^4+(-I*Psi^6*a[1]+2*Psi^5*a[1]+((2*I)*a[1]+I*b[1])*Psi^4+(-4*a[1]+4*b[1])*Psi^3+(-I*a[1]-(2*I)*b[1])*Psi^2-2*b[1]*Psi+I*b[1])*cosh(d)^2-Psi^2*(I*Psi^2*a[1]+(-2*a[1]+2*b[1])*Psi-I*b[1]))*(Psi+1)*abs(Psi(t))^2*(Psi-1))/(abs(psi(t))*cosh(d)^4):
Eqs := {seq(EQ[i], i = 0 .. 6)};
       {EQ[0], EQ[1], EQ[2], EQ[3], EQ[4], EQ[5], EQ[6]}

@dharr 

I vote up for all the impressive work you did here.

@dharr 
(I reply from my personal account)

I'm not very familiar with the notion of the square root of a matrix, I use it in some statistics algorithms and it's almost all.
As I often manipulate big correlation matrices (of floats), I looked for algorithms that would be cheaper, in time and/or memory, than the eigenvalues decomposition.
A lot of them are based on a generalization to matrices of the Newton method used for iteratively assessing the square root of a number.
I've coded two of these iterative algorithms: they are almost as faster than the method based on eigenvector decomposition, but they used about 1% of the memory size this latter requires (3 Gb for a 100x100 dense correlation matrix).

But the main problem is that these methods do not always give the same square root, even when the matrix is positive-definite (the square root matrix must then be unique).
Starting from a matrix S (symmetric definite positive) and computng the square root matrix Q of C=S^2, I can see that the eigenvalues decomposition finds Q=S, and so does one Newton-like method, but the other one finds a matrix Q extremely different from S, although Q^2=C.
I can't figure out where this nonsense come from.
Did you know about these Newton-like methods?
If it is so do you have any idea about this issue?

If you think this issue deserves a new thread I can ask a specific question, let me know

Download MatrixSquareRoot.mw

@AHSAN 

Why do you use the shooting method to solve your problem? I don't know, I'm not familiar enough with this method to give you a reliable answer.

@Carl Love ?

I read your first response and I don't know if you are wrong, maybe it's me.
But, to be honest,I didn't care about the confusing code provided by the OP.  What I coded are the equations the OP presented before the code, the ones he says he wants to solve and for which bn is obviously a parameter.

I updated my answer to do something which, I think, corresponds to what the OP wants.
Basically I believe the OP has 2 ode systems parameterized by bn and he wants to plot its solution wrt bn.
For instance, as z is a function of x, the plot could be either a set of curves (x, z(bn)(x)) for different values of bn, either a surface (x, nb, z(bn)(x))

Waiting for his correction and complements.

@MaPal93 
 

Answer to your first question.
As the output exceeds 1000000 you can't display it if you use collect in a raw way.
So here is an indirect way to collect any components of the solution wrt to any random variables. 
I am not sure that the final expressions present any interest given their complexity, but that is my opinion and the decision is yours.
290423_Chi_version1_mmcdara.mw

Answer to your last question.
Here is a notional example
Last_question.mw

@Anthrazit 

I guess here is the reason why the two plots do not have the same height:

For some reason I don't know geometry:-draw displays CONSTRAINED scaled object, and this property is seems to be shared by all the plots that are displayed with a drawn object... unless you explicitly declare SCALING(UNCONSTRAINED).

Download The_reason_is.mw

BTW: you're right, panning upwards or downwards displays the correct figure.

@Anthrazit 

You wrote "When panning the graphics on the right side the plot is clipped."
Here is WHAT I SEE  when I open your file with Maple 2015:

And this is WHAT I GET when I execute the worksheet (with Maple 2015) and WHAT YOU SHOULD SEE when you open the file I sent you with Maple 20xx (before executing it!)

This is why I said there was no problem with Maple 2015.

No problem with Maple 2015.2 (even if I don't understand why the plots appear befor the command lines and not after)
plotpoint2_2015.mw

@AHSAN 

I've just updated  my first answer

@Jamie128 

LOL

I guess you are aware that 

combine(convert(A1, trig), trig)
                            cos(theta) 
                         - ------------
                           2 sin(theta)

?

Have you try using the DeepLearning package?
https://www.maplesoft.com/support/help/maple/view.aspx?path=DeepLearning

@shreen 

C1 is the output of 

C1 := C_Solution[1]:

Thus 

C2 := C_Solution[2]:

and so on.
Can't you do it yourself?

Have you seen the size of the expression which represents C1???
What do you think that tits display will reveal to you? Look at this:
calc_slip_two_cylinders_electo_mmcdara.mw

@acer 

Thank you acer for having corrected this.

[ A followup: 4eqns.mw ]